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While  the  new  theoretical  treatment  given  here  to  the  finite- 
element  method  resulted  in  discovering  an  optimum  finite-element 
procedure  which  is  actually  a special  case  of  a finite-difference  method 
which  was  reported  by  Crandall,  it  lays  the  framework  for  a similar 
approach  which  can  be  used  to  give  a truly  useful  numerical  tool:  a 
finite-element  procedure  for  parabolic  problems  which  results  in  a 
higher  order  of  accuracy  in  the  mean  square  sense  over  the  entire 
solution  domain. 
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of  the  Air  Force  Institute  of  Technology,  for  his  advice  on  special 
occasions,  and  to  Dr.  W.  Kessler,  of  the  Air  Force  Materials  Laboratory, 
for  sponsoring  this  research  project.  I wish  also  to  express  my 
gratitude  to  my  wife,  Susann,  who  helped  in  preparation  of  the  final 
draft  and  gave  invaluable  moral  support  all  during  the  project,  and 
to  Sharon  Flores,  whose  expert  typing  skills  added  much  to  this  thesis. 
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Abstract 

The  transient  heat  conduction  equation,  with  Dirichlet  and  Neumann 
boundary  conditions,  is  solved  by  the  methods  of  finite-differences 
and  finite-elements,  and  the  numerical  solutions  are  investigated  with 
respect  to  accuracy  and  stability.  A general  six-point  finite- 
difference  expression  is  used  for  which  there  exists  a high  order 
accurate  modification.  The  version  of  the  finite-element  method  used 
was  based  on  a variational  principle  which  is  stationary  in  time;  the 
temporal  behavior  of  the  differential  equation  is  treated  with  a 
finite-difference  approximation.  This  method  is  shown  to  be  equivalent 
to  the  method  of  Galerkin.  Several  methods  for  treating  accuracy  and 
convergence  problems  which  result  from  a discontinuity  in  the  initial 
condition  are  investigated.  Hie  Crank-Nicolson  method  is  a special 
case  of  both  the  finite-difference  and  finite-element  methods.  The 
finite-difference  version  of  the  Crank-Nicolson  method  is  shown  to  be 
more  accurate  than  the  finite-element  version,  especially  when  a dis- 
continuity exists  between  the  initial  condition  and  the  boundary 
conditions.  The  high  order  accurate  schemes  for  both  finite-differences 
and  finite-elements  are  shown  to  be  equivalent  for  the  case  of  linear 
elements  derived  from  a stationary  variational  principle.  Some  of  the 
results  suggest  the  possiblity  of  finding  a finite-element  scheme  which 
is  highly  accurate  in  a mean  square  sense  over  the  entire  soltuion 
domain. 
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AN  INVESTIGATION  OF  THE  NUMERICAL  METHODS 
OF  FINITE  DIFFERENCES  AND  FINITE  ELEMENTS 
FOR  DIGITAL  COMPUTER  SOLUTION  OF 
TOE  TRANSIENT  HEAT  CONDUCTION  (DIFFUSION) 
EQUATION  USING  OPTIMUM  IMPLICIT  FORMULATIONS 

I.  Introduction 


Background 

Many  problems  arise  in  engineering  practice  which  involve  partial 
differential  equations.  Exact  analytical  methods,  such  as  the  classical 
separation  of  variables  technique,  often  cannot  be  used  to  solve  these 
equations.  The  purpose,  then,  of  numerical  techniques  is  to  obtain 
reasonably  accurate  results  for  those  problems  for  which  exact  methods 
cannot  be  used. 

The  most  well  known  numerical  technique  is  the  finite-difference 
method.  This  method  approximates  the  derivatives  which  occur  in  partial 
differential  equations,  reducing  them  to  a set  of  algebraic  equations 
which  are  then  solved  on  a digital  computer. 

During  the  1960's,  there  was  widespread  development  of  a new 
method,  the  finite-element  method,  for  obtaining  numerical  solutions 
to  differential  equations.  Most  of  these  developments  were  in  the  field 
of  solid-body  mechanics,  but  more  recently  there  have  been  applications 
of  the  finite-element  method  in  the  area  of  heat  transfer.  The  finite- 
element  method  converts  the  original  partial  differential  equation 
into  a variational  integral  which  must  be  minimized;  the  solution  of 
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The  original  partial  differential  equation  will  minimize  this  integral. 

As  with  the  finite-difference  method,  application  of  the  finite-element 
method  results  in  a set  of  algebraic  equations  which  can  be  solved 
on  a digital  computer. 

One  application  of  these  numerical  methods  lies  in  the  area  of 
heat  transfer.  It  is  part  of  the  responsibility  of  the  Air  Force 
Materials  Laboratory  to  analyze  the  thermal  response  and  ablation 
characteristics  of  rocket  nozzles  and  re-entry  vehicles.  The  transient 
heat  conduction  equation,  a boundary  value  problem,  is  encountered  in 
this  analysis.  Because  of  the  irregular  geometry  involved,  approximate 
nunerical  techniques,  such  as  finite-differences  and  finite-elements, 
must  be  used. 

One  of  the  more  accurate  techniques  of  approximating  the  transient 
heat  transfer  equation  is  the  Crank-Nicolson  method  (Ref  1).  This 
technique  was  developed  originally  for  use  with  the  finite-difference 
method.  There  is  a method  which  is  related  to  the  Crank-Nicolson  scheme, 
the  Crandall  Method  (Ref  2),  which  minimizes  the  truncation  error  in 
the  finite-difference  method.  There  is  also  a version  of  the  Crank- 
Nicolson  scheme  which  was  developed  for  use  with  the  finite-elements 
method.  Obviously,  if  the  error  in  the  finite-element  method  could  be 
minimized  in  a manner  similar  to  the  Crandall  method,  then  the  resulting 
solution  would  be  more  accurate.  Finding  this  optimum  method  was  the 
ultimate  goal  of  this  research  project.  The  increased  accuracy  could 
be  very  beneficial  to  heat  transfer  engineers,  provided  it  does  not  come 
at  the  expense  of  some  more  important  factor,  such  as  computer  run  time 
or  stability. 
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Problem 

The  objective  of  this  project  was  to  solve  a one-dimensional 
version  of  the  transient  heat  conduction  equation,  with  a known 
analytical  solution,  using  modifications  of  the  Crank-Nicolson  versions 
of  the  finite-difference  method  and  the  finite-element  method.  The 
objective  was  to  obtain  the  most  accurate  and  stable  solution  to  the 
differential  equation.  Accuracy  is  defined  here  as  the  deviation  of 
the  approximate  solution  from  the  exact  analytical  solution.  Stability, 
on  the  other  hand,  refers  to  unwanted,  numerically  induced  oscillations 
which  may  affect  the  solution.  An  unstable  condition  results  when 
these  oscillations  grow  without  bound. 

Scope 

Initially,  the  analysis  was  only  to  include  the  following:  (1)  a 
comparison  of  the  Crank-Nicolson  version  of  the  finite-difference  method 
with  the  Crank-Nicolson  version  of  the  finite-element  method  with 
respect  to  accuracy  and  stability,  and  (2)  an  attempt  to  reduce  the 
error  in  the  finite-element  method  by  experimentally  varying  a parameter. 
The  analysis  was  not  initially  planned  to  include  a theoretical  treatment 
of  the  optimization  process.  However,  during  the  investigation,  the 
scope  was  broadened  to  include  the  explicit  (Euler)  and  fully  implicit 
methods,  and  also  a theoretical  treatment  of  the  optimization  process 
for  finite- e 1 ements . 

Assumptions 

The  fundamental  assumptions  which  will  hold  throughout  the  text 
are:(l)  that  the  physical  properties  of  the  material  of  interest  are 
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constant  with  respect  to  time  and  spatial  variables,  (2)  that  the 
analysis  in  one  dimension  is  representative  of  reality  (for  example, 
an  insulated  thin  rod  or  a plane  wall  can  be  adequately  represented 
with  one  spatial  dimension  under  certain  conditions),  (3)  that  constant 
mesh  spacing.  Ax  , is  adequate  for  the  user's  needs,  and  finally, 

(4)  that  no  heat  generation  takes  place  within  the  material  of  interest. 

It  is  obvious  that  the  restrictions  are  severe.  One  of  the 
biggest  difficulties  encountered  when  using  the  finite-difference 
method  occurs  when  either  the  arrangement  of  the  nodes  or  the  geometry 
is  irregular.  Some  of  these  problems  are  simplified  by  the  finite- 
element  method.  Unfortunately,  these  advantages  may  be  negated  by  the 
restriction  to  constant  material  properties  and  constant  mesh  spacing. 
But  the  gains  in  accuracy  are  sure  to  find  applications  even  under 
these  restrictions.  Also,  it  should  be  noted  that  this  analysis  does 
not  necessarily  hold  for  two  or  three  spatial  dimensions,  although  it 
is  the  opinion  of  the  author  that  the  method  can  be  extended. 

The  restriction  on  heat  generation  here  is  merely  for  convenience, 
and  the  results  of  this  paper  would  not  be  affected  by  the  inclusion 
of  a heat  generation  term  so  long  as  this  term  is  constant  with  respect 
to  the  time  and  spatial  variables. 
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II.  Theory 

The  Physical  Problem 

General.  The  linear  partial  differential  equation  of  second 
order  in  u(w,v) 


Au  + Bu  + Cu  + Du  + Eu  + Fu  = G (1) 
ww  wv  w w v 

where  A,B,...,G  are  constants,  or  functions  of  w and  v only,  and 
where  the  subscript  indicates  differentiation  with  respect  to  the 
subscript,  is  classified  as  hyperbolic,  parabolic,  or  elliptic  type  in 
a domain  of  the  (w,v)  plane  depending  on  whether  the  values  of  the 
discriminant  B2-4AC  are  positive,  zero,  or  negative,  respectively, 
throughout  the  domain. 

The  heat  conduction  equation  is  parabolic  in  form  and  is  no  more 
than  a mathematical  statement  of  the  first  law  of  thermodynamics  and 
Fourier's  Law  of  heat  transfer.  The  first  law  of  thermodynamics 
states  that  the  overall  change  in  the  internal  energy  of  a system  is 
equal  to  the  heat  which  flows  into  the  system  plus  the  heat  generated 
within  the  system  less  the  heat  which  flows  out  of  the  system.  If 
the  derivative  of  each  of  these  quantities  is  taken  with  respect  to 
time,  the  result  is  the  following  energy  balance: 


0 


STOR 


0.  + 0 

in  gen 


(2) 
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where  the  dot  indicates  differentiation  with  respect  to  time. 

The  energy  terms  in  Equation  (2)  are  replaced  by  the  appropriate 
rate  equations.  Some  of  the  more  common  rate  equations  are  for 
conduction,  convection,  radiation,  heat  storage,  and  heat  generation. 
Conduction  is  defined  as  heat  transfer  through  matter  without  net  movement 
of  the  material,  such  as  in  a solid.  The  rate  equation  is  given  by 
Fourier's  Law: 


where 


-kA§I 

dX 


■ rate  of  heat  flow  in  the  x direction,  Btu/hr 

* coefficient  of  thermal  conductivity,  Btu/hr-ft-°F 

■ area  normal  to  the  x direction  through  which  heat  flows,  H2 
» temperature,  #F 

= space  variable,  ft 


Heat  storage  is  described  mathematically  by 


^STOR  * pVcp  3t 


where 


rate  of  heat  storage,  Btu/hr 
density,  lbm/ft3 


V * volume,  ft3 

Cp  ■ specific  heat,  Btu/lbm-#F 

t * time,  hr 

The  expressions  for  convection,  radiation,  and  heat  generation  will 
not  be  needed  in  this  report,  but  may  be  found  in  Myers  (Ref  3:2-4). 


Figure  1.  A Unit  of  Volume  From  a Wall  with  Large 
Dimensions  in  the  y and  2 Directions. 


From  Figure  1,  the  energy  balance  equation  for  the  given  plane 
wall  system  can  be  seen  to  be 


»ca|I 

P 3t 


- u 57 

9 x+  sx  9x  x 
fix 


(5) 
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where  it  is  assumed  that  no  energy  generation  takes  place  within  the 
unit  of  volume  shown,  and  A represents  the  unit  of  area  shown.  In 
the  limit  as  Ax  approaches  zero,  the  equation  becomes  the  standard 
heat  conduction  equation  which  is  parabolic  in  type: 


* 3T  3 ri.  * 

pcpA  at  = T7  tkA 


3x 


3T, 

3xJ 


(6) 


If  the  material  properties  p , Cp  and  k are  constant  with  respect 
to  x , and  a constant  cross-sectional  area  is  used,  the  equation 
becomes 


II  JL. 

8*  ' “W 


(7) 


This  is  the  form  of  the  heat  equation  which  will  be  analyzed  in 
this  paper.  The  problem  is  not  completely  specified  until  the  heat 
equation  is  associated  with  boundary  conditions  and  initial  conditions. 

Two  types  of  boundary  conditions  will  be  investigated,  Neumann  and 
Dirichlet.  If  the  function  itself  is  specified  on  the  boundary,  it  is 
termed  a Dirichlet  condition.  If  the  derivative  of  the  function  is 
specified  along  the  boundary,  it  is  termed  a Neumann  condition.  If 
the  boundary  condition  contains  values  of  the  function  and  its  derivative, 
it  is  referred  to  as  a mixed  boundary  condition. 

Primary  Problem.  The  primary  problem  considered  is  that  of  a 
parallel  sided  plane  wall,  infinite  in  all  directions  normal  to  the 
direction  of  heat  flow,  which  is  heated  until  a steady  state  temperature 
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of  1000°  F is  reached  at  all  points  within  the  wall.  The  surfaces 
of  the  wall  are  then  cooled  instantaneously  to  0°  F and  maintained  at 
that  temperature.  This  problem  was  chosen  for  its  simplicity  and  also 
because  it  has  been  analyzed  by  the  Crank-Nicolson  and  Crandall  finite- 
difference  methods  (Ref  4:325).  In  order  to  provide  more  generality  to 
the  problem,  it  will  be  analyzed  in  normalized  form.  If  Tg  is  the 
specified  boundary  temperature  and  T^  is  the  initial  temperature, 
then  the  problem  can  be  written 


3T 

at 


a2T 
a 

“ 3x2 


(8) 


T(x,t)  * T.,  t = 0 (9) 

T(0,t)  = T(L,t)  =■  TB,  t > 0 (10) 


where 

am  * k/pc^  = the  thermal  diffusivity  of  the  wall  material 

L ■ the  thickness  of  the  wall 

Ta  * 1000°  F 

Tb  - 0“  F 
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Let  u ■ 


(T-TJ/T  , x * x/xn  , and  0 - t/t  where  T 

BOO  O O 

x , and  t are  left  unspecified  for  the  tine  being.  These  dimension 
o o 

less  variables  are  then  substituted  into  Equation  (8)  which  is  then 
multiplied  by  tQ/xQ2  to  yield 


3u 

30 


t 


x 


o 


o 32u 

2 X 2 


(11) 


Since 

t 

o 


xq  and  tQ  are  unspecified,  the  choice  of  xq  = L and 

L2/a  yields  the  desired  normalized  partial  differential  equation: 

01 


3u 

30 


If  T is  chosen  as  T. 
o i 

after  substitution 


3fu 

3x2 


(12) 


the  boundary  and  initial  conditions  become. 


u(x,0)  = 1 ,0  = 0 (13) 

* 

u(O,0)  « u(l,0)  *0  , 0 > 0 (14) 

A schematic  diagram  of  the  normalized  problem  is  shown  in  Figure  2. 
The  most  striking  feature  of  the  problem  is  the  discontinuity  between 
the  initial  condition,  u(x,0)  = 1 , and  the  boundary  conditions, 

u(O,0)  ■ u(l,0)  = 0 . This  discontinuity  affects  the  solution  most 

noticeably  at  points  near  the  boundary  and  at  early  times,  and  it  is 
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Figure  2.  A Schematic  Diagram  of 
the  Primary  Problem. 

the  source  of  problems  with  convergence  (to  be  discussed  later)  of 
both  the  finite-difference  and  finite-element  higher  order  accurate 
schemes.  Several  attempts  to  reduce  the  effect  of  this  discontinuity 
on  the  numerical  solution  of  the  problem  will  be  discussed. 

The  exact  solution  of  this  problem  can  be  obtained  by  a number 
of  techniques.  The  method  of  separation  of  variables  was  used  here, 
but  the  Laplace  transform  method  yields  the  same  form  of  solution. 

The  problem  is  solved  in  Appendix  A yielding  an  infinite  series  as  the 
solution: 

u(x,0)  = l sin((2n-l)n  x)  e”^2n_1^  9,  6 > 0 

n=l 


The  solution  is  shown  in  graphical  form  in  Figure  3. 


(15) 
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TEMP  (NORMALIZED) 

J).00  0.17  0.33  0.50  0.67  0.83  1.00  1.17 

rj i i i i i i 


Verification  of  the  Primary  Problem.  The  coefficients, 

a ■ 4/(2n-l)n)  , are  bounded  for  all  n since  a -*•  0 as 
n n 

n -*■  “ ; therefore,  la  I < M , where  M is  some  constant.  Thus, 

n ' 

whenever  0 2 0q  , where  0q  is  any  positive  constant 

/ / 

|an  sin((2n-l)n  x)  exp(-((2n-l)*)20|  < M exp(- ( (2n-l)n)20Q)  (16) 

Since  the  ratio  test  ensures  the  convergence  of  the  infinite  series 

with  constant  terms  M exp(-((2n-l)ir)20o)  , the  Weierstrass  M-test 

provides  a sufficient  condition  for  the  uniform  convergence  of  Equation 

(15)  when  0 f 3T  < 1 , 0 J 0q  > 0 (Ref  5:28).  Since  the  terms  of 

the  series  in  Equation  (15)  are  continuous  functions,  the  series 

_ 

converges  to  the  continuous  function  u(x,0)  for  0 > 0 , since 

0q  is  any  positive  constant. 

Similarly,  the  series  with  terms  exp(-((2n-l)n)20o)  , or 

(2n-l)w  exp(-((2n-l)it)20o)  , also  converges.  Therefore,  the  series 

in  Equation  (15)  can  be  differentiated  once  with  respect  to  t and 
twice  with  respect  to  x , when  t > 0 , since,  by  the  Weierstrass 

M-test,  the  series  of  derivatives  converges  uniformly.  Further,  since 
the  terms  of  Equation  (15)  satisfy  Equation  (12),  by  superposition, 
the  sum  of  the  series,  u(x,0)  , also  satisfies  Equation  (12). 

To  establish  the  convergence  of  the  series  in  Equation  (15)  to 
the  initial  condition  at  0*0  , use  will  be  made  of  Abel’s  test 

(Ref  5:228).  First,  f(T)  *s  defined  as  the  odd  periodic  extension  of 
the  initial  condition  with  a period  of  2.  Further,  for  each  xeR  , 
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the  series  with  terms  an  sin((2n-l) w x)  converges  to  f(x)  . At 
points  where  a discontinuity  exists,  f(x)  is  defined  as 

f(J)  . f(x  . 0)  ♦ ttx  r 0)  m 

where  f(x  + 0)  and  f(3T  - 0}  are  the  right  and  left  hand  limits  of 
f at  x . According  to  Abel's  test,  the  series  which  results  after 
multiplying  each  term  of  a convergent  series  by  the  corresponding 
elements  of  a sequence  of  functions  of  0 , which  is  bounded  from 

above,  as  is  exp(-((2n-l)ff)29)  , is  uniformly  convergent  for  all 

0 > 0 . 

This  completes  the  formal  verification  of  the  solution  to  the 
primary  problem. 

Error  in  Truncation  of  the  Infinite  Series.  The  computer, 
although  fast,  is  incapable  of  summing  an  infinite  series.  However, 
it  would  be  of  value  to  estimate  the  error  which  results  when  the 
infinite  series  solution  is  truncated  after  N terms.  The  following 
theorem,  given  by  Thomas  (Ref  6:660),  gives  the  conditions  necessary  to 
estimate  the  error  incurred  by  truncating  an  infinite  series. 

Theorem.  If  the  series 

(18) 


! 


I a sin((2n-l)*  x)  exp(-((2n-l)ir)20) 

n*l 

is  strictly  alternating,  and  if  the  n^1  term  tends  to  zero  as 
n -*■  • , then  each  term  is  numerically  less  than  or  equal  to  its 
predecessor.  The  proof  of  this  theorem  is  given  by  Thomas.  One 
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further  detail  is  required  to  get  an  error  estimate.  The  terms  of 
the  series  (18)  must  be  regrouped  into  a strictly  alternating  series. 
This  can  be  done  if,  and  only  if,  (18)  is  uniformly  convergent,  which 
has  previously  been  shown  to  be  the  case.  Inspection  of  (18)  indicates 
the  second  condition  above  is  met,  also.  The  regrouping  is  done  by  the 
computer  by  testing  the  sign  of  each  term  of  the  series  sequentially 
and  adding  the  current  term  to  the  previous  term  if  the  sign  remains 
unchanged,  or  by  beginning  a new  term  if  the  sign  changes.  The  error 
estimate  is  given  by  the  last  complete  term  of  the  new  strictly 
alternating  series.  The  series  terminates  when  the  computer  reaches 
the  limits  of  storage  capacity  prior  to  an  underflow  condition,  that  is 
when  a number  is  so  small  that  the  computer  cannot  represent  it  in  a 
memory  word. 

Of  course,  the  attention  to  detail  here  is  somewhat  academic, 
since  the  sine  and  exponential  functions  used  in  the  analytic  solution 
are  themselves  approximated  by  built  in  computer  functions  which  make 
use  of  truncated  infinite  series. 

Secondary  Problem.  The  secondary  problem  considered  is  not  given 
as  complete  a treatment  as  the  primary  problem.  The  only  difference 
between  the  two,  in  the  unnormalized  form,  is  that  the  right  boundary 
temperature  is  not  fixed  for  all  t > 0 , but  rather  that  boundary  is 

insulated.  Therefore,  the  temperature  gradient  at  the  right  boundary 
is  zero.  Normalization  of  the  secondary  problem  is  similar  to  that  of 
the  primary  problem  and  yields 


3u 

36 


(19) 
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The  normalized  problem  is  shown  schematically  in  Figure  4. 

Again,  there  is  a discontinuity  between  the  initial  condition  and  the 
left  boundary  condition.  The  exact  analytical  solution  of  the  secondary 
problem  is  given  in  Appendix  B.  TT\e  method  of  separation  of  variables 
was  again  used  with  the  following  result: 


U(X,0)  * \ (2n-i')V  »in((2n-l)J)  e‘ t(2n_1)Il2e  (23) 

Verification  of  the  solution  is  exactly  analogous  to  that  of  the 
primary  problem  and  will  not  be  repeated  here.  The  solution  appears 
in  graphical  form  in  Figure  5. 

Finite-Difference  Formulation 

General  Expression.  The  finite-difference  method  makes  use  of 
differences  to  approximate  derivatives.  For  example 


du  „ % ~ ua 
dx  ' \ - x. 


(24) 


could  be  used  to  approximate  the  derivative  of  u with  respect  to 
x taken  at  point  b,  point  a,  or  halfway  between  a and  b.  If  it 
represents  the  derivative  at  a in  Figure  6,  it  is  a forward  difference 
expression.  If  the  approximation  represents  the  derivative  at  point  b, 
it  is  a backward  difference  expression,  and  if  it  approximates  the 
derivative  at  the  point  halfway  between  a and  b,  it  is  a central 
difference  expression.  The  central  difference  approximation  is  the 
most  accurate. 
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Figure  6.  Finite-Difference  Approximation 
of  a Derivative. 


These  expressions  can  be  used  to  replace  the  derivatives  in  the 
normalized  heat  equation.  If  the  heat  equation  is  forward  differenced 
in  time  and  central  differenced  in  space,  the  resulting  expression  is 


ui,k+l  ~ Ui,k 
A0 


u.  , . - 2u.  , + u.  . . 
l-l  ,k i,k  i+l,k 

AX2 


(25) 


where  the  subscript  i refers  to  a nodal  point  on  the  space  axis, 
and  the  subscript  k refers  to  a point  on  the  time  axis  in  Figure  7. 
The  bar  over  the  x has  been  dropped  for  simplicity.  This 
expression  is  called  explicit  sin^e  the  temperature  at  the  new  time, 
ui  k+1  ' *s  solve<*  explicitly  in  terms  of  temperatures  on  the 

previous  time  level,  k 
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Figure  7.  Space-Time  Grid  for 
the  Heat  Equation. 

If  the  heat  equation  is  backwards  differenced  in  time  and  central 
differenced  in  space,  the  expression  is 


ui,k»l  ~ uj.k 
A0 


Ui-l,k+l  " 2uj ,k+l  * Uj+l,k+l 
Ax2 


This  expression  is  termed  implicit,  since  it  requires  the  solution  of 
a set  of  simultaneous  equations  at  each  new  time  level,  shown  here 
schematically  as  k+1  . 

In  1947,  Crank  and  Nicolson  (Ref  1)  proposed  a scheme  which  can 
be  considered  as  a central  difference  expression  in  time  at  time  level 
k + j , and  it  makes  use  of  central  differences  in  space.  The  result 
is  a more  accurate  difference  equation: 


(26) 
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ui,k+l  ~ ui,k 
A0 


1 ui-l,k+l  ~ 2uj,k+l  * ui+l,k+l 

2 Ax2 


(27) 


1 ui-l,k  " 2ui ,k  + Ui+l,k 
+ 2 Ax2 


This  is  also  an  implicit  expression,  as  it  requires  the  solution  of  a 
set  of  simultaneous  equations.  Note  the  equal  weight  on  the  new  time 
level,  k+1  , and  the  old  time  level,  k , on  the  right  hand  side 

of  the  equation. 

In  1955,  Crandall  (Ref  2:318)  proposed  a generalized  expression 
replacing  the  equal  weights  shown  in  Equation  (27)  with  variable  ones 
as  shown  in  Figure  7.  The  resulting  expression  was 


Ui,k+1  " ui,k 
A0 


ui-l ,k+l  " 2ui,k+l  * ui+l,k+l 
Ax2 


♦ (l-o) 


Ui-l,k  " 2ui,k  * ui+l,k 
AX2 


(28) 


Crandall  observed  that  for  a certain  value  of  the  parameter  a , 
the  expression  becomes  very  accurate.  For  this  optimum  value  of  a , 
the  error  incurred  in  approximating  the  time  derivative  very  nearly 
equals  the  error  in  the  approximation  of  the  space  derivative.  The 
optimum  value  of  a , shown  in  Figure  8,  is  given  by 
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THEORETICAL  OPTIMUM  ALPHA 
VALUES  FOR  FINITE  DIFFERENCES 
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FOURIER  MODULUS  (P) 


2.00 


Figure  8,  Theoretical  Optimum  Alpha  Values  for  Finite  Differences. 
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a ■ 


(29) 


i.  (i L) 

2 U 6pJ 


where  p * A0/Ax2  . The  dimensionless  parameter,  p , is  known  in 
the  literature  as  the  Fourier  modulus  (Ref  19:938). 

Application  of  the  General  Expression.  To  apply  the  general 
expression.  Equation  (28),  to  the  primary  problem,  the  space  axis  in 
Figure  7 is  first  divided  into  N-l  intervals.  The  value  of  the  first 
node  is  governed  by  the  boundary  condition  at  x=0  ; whereas,  the 

value  of  the  Nth  nodal  point  is  governed  by  the  boundary  condition 
at  x»l  . As  the  subscript  i is  varied  over  the  range  of  nodal 
values,  the  following  set  of  equations  results,  written  here  in  matrix 
form  for  the  case  N-6  : 


0 

l+2pa 

-pa 

0 

0 

0 

-p  a 

l+2pa 

-pa 

0 

0 

0 

-pa 

l+2pa 

-pa 

0 

0 

0 

-pa 

l*2pa 
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l-2p(l-a) 

P(l-a) 

0 

0 

p(l-a) 

l-2p(l-a) 

P(l-a) 

0 

0 

P(l-a) 

l-2p(l-a) 

Pd-a) 

0 

0 

p(l-a) 

l-2p(l-a) 

or  in  vector  notation 


A u a B u (31) 


Equation  (30)  has  four  equations,  but  six  unknowns.  This 
difficulty  is  avoided,  however,  since  the  end  nodes  are  not  needed  and 
can  be  discarded.  The  result  is 


■ J-V  ’ -'V 


l-2p(l-a) 

P(l-a) 

0 

Pd-c0 

l-2p(l-a) 

p(l-a) 

0 

p(l-a) 

l-2p(l-a) 

0 

0 

p(l-a) 

0 

U2 

0 

U3 

Pd-a) 

U4 

l-2p(l-a) 

us 

(32) 


Note  that,  because  of  symmetry,  only  two  columns  of  information  will 
need  to  actually  be  stored  in  computer  memory  by  a production  computer 
code  for  each  of  the  coefficient  matrices.  Additionally,  the  matrices 
are  tridiagonal;  thus,  a matrix  factorization  scheme,  the  Thomas 
method,  may  be  used  to  solve  the  matrix  equation.  The  Thomas  method 
(Ref  7:46-48)  has  the  advantage  of  requiring  fewer  algebraic  operatons 
than  the  more  familiar  Gauss  reduction  method.  There  are  also  iterative 
methods  for  solving  these  systems  of  equations. 

To  apply  Equation  (29)  to  the  secondary  problem,  a fictitious  point 
must  be  created  where  node  N+l  would  be  if  it  existed  (Ref  8:35). 
Applying  Equation  (28)  to  node  N yields 


* (1*2p“>  “N.k.l  - poVl,k*l 


■ Pt'-^Vl.k  * 


(33) 

* P(l-«>Vl,k 
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i-2P<i-“>Yk 


Additionally,  the  derivative  boundary  condition  at  x»l  must  be 
approximated . Thus 


0 , 6 > 0 


becomes 


^Wl.k  " uN-l,k 


where  a central  difference  approximation  has  been  used.  Solving 
Equation  (35)  for  u^j  ^ and  substituting  back  into  Equation  (33) 
yields 


>Vl,Kl  * (1*2>X,)Vk.l  ■ 2p(‘-<‘)UN-l.k  * (l-2PCl-o))uN>k 


Again  using  N=6  , the  matrix  equation  now  becomes 


a-  • 


I- 


l*2pa 

-pa 

0 

0 

0 

-Pa 

l+2pa 

-pa 

0 

0 

0 

-pa 

l*2pa 

-pa 

0 

0 

0 

-pa 

l*2pa 

-pa 

0 

0 

0 

- 2 pa 

l+2pa 

U2 

U3 

U4 

U5 

1 

c 

0 

1  

k+1 


l-2pd-o)  p(l-a) 

p(l-a)  l-2p(l-a) 


p(l-a) 
p(l-a)  l-2p(l-o) 


Pd-a) 
P(l-a)  l-2p(l-o) 


p(l-a) 
2p(l-a)  l-2p(l-a) 


"*  *" 

u2- 

U3 

U4 

U5 

1 

vO 

3 

1 

Linear  Finite  Element  Formulation 


Background . The  finite-element  method  is  fundamentally  different 
in  nature  from  the  finite-difference  method.  The  finite-difference 
method  is  designed  to  approximate  the  solution  of  a differential 
equation  at  a number  of  specified  nodal  points;  this  method  gives  no 

27 


(37) 
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information  about  the  solution  between  these  nodal  points.  The 
finite-element  method  does,  however,  assume  a general  form  for  the 
solution  between  these  nodal  points.  The  method  will  be  introduced  by 
way  of  the  following  simple  illustration. 

The  differential  equation 


♦ 1)  u * f 


u(0)  ■ u(L)  * 0 


bears  a special  relationship  to  the  functional 


1(0)  - j [DO.Q)  - 2(f,0)] 


where  the  notation  (gj.gj)  represents  the  inner  product  of  the 
two  functions  g^  and  g2  , or 


(gj'Sj)  * I gxM  g2(x)  dx 


and  0(x)  ■ u(x)  ♦ £v(x)  . The  function  0 is  called  a trial 
function,  and  the  only  restriction  on  the  function  v , called  a 
variation  of  u , is  that  v must  satisfy  all  Dirichlet  boundary 
conditions,  those  written  in  terms  of  the  function  itself.  Thus, 
v must  have  the  properties 


j 


j 

v(0)  - v(L)  « 0 (42) 

such  that 

f 

/ 

0(0)  * u(0)  = Q(L)  = u(L)  = 0 (43) 

Except  for  this  restriction,  v(x)  is  left  completely  arbitrary.  Now 

the  relationship  between  Equation  (38)  and  Equation  (40)  can  be 

demonstrated.  The  functional  1(0)  has  a minimum  value  at  the  point 

where  0 = u . That  is,  at  the  point  where  £ = 0 . Therefore,  if 

the  minimum  value  of  the  integral  1(0)  can  be  found,  the  value  of  0 

which  minimizes  I is  equal  to  the  solution,  u , of  Equation  (38). 

In  the  finite-element  method,  the  domain  of  the  solution  is 

divided  into  N elements.  Within  each  of  these  elements,  the  temperature 

distribution  is  assumed  to  have  a certain  shape,  for  example,  linear, 

as  in  Figure  9.  A discussion  of  the  notation  and  the  meaning  of 

certain  terms  would  be  helpful  at  this  point. 

The  temperature  distribution,  u^  , where  the  superscript 

indicates  the  number  of  elements  in  the  solution  domain,  is  a member 

N N 

of  the  finite-dimensional  space  of  functions,  S ; the  space,  S , 
is  a subspace  of  the  infinite  dimensional  Hilbert  space,  H*  , which, 
in  addition  to  the  general  defining  properties  of  a Hilbert  space,  is 
defined  as  the  space  of  functions  which  are  square  integrable  and  whose 
first  derivatives  also  are  square  integrable.  Further,  the  subscript 

\ 
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Figure  9.  Finite-Element  Arrangement 
for  Solution  of  Eq.  (38), 

indicates  that  only  the  essential  boundary  conditions,  defined  below, 
of  Equation  (39)  must  be  satisfied  (Ref  9:5-11).  In  this  particular 
case,  both  conditions  are  essential. 

An  essential  boundary  condition  is  one  whose  statement  is  in  terms 
of  derivatives  below  order  s when  the  differential  equation  is  of 
order  2s  . The  reason  for  this  distinction  between  essential 
conditions  and  natural  conditions,  which  are  defined  to  be  those  written 
in  terms  of  derivatives  or  order  s and  higher,  lies  at  the  heart 
of  the  method.  The  inhomogeneous  term,  f , in  Equation  (38)  will  be 
used  to  illustrate  convergence  in  the  mean  in  hopes  of  clarifying  the 
distinction  made  between  essential  and  natural  boundary  conditions. 


!i 
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The  eigenfunction  expansion  of  f is  made  in  terms  of 
functions  which  necessarily  satisfy  the  boundary  conditions  on  the 
problem.  However,  the  function  f itself  does  not  necessarily 
satisfy  those  conditions.  The  inifinte  series  expansion  of  f in 
V,  terns  of  the  eigenfunctions  is  said  to  converge  in  the  mean  square 
\*gnse  to  f if 


where 


/ [f(x)  - Sm(x)]2  dx  = 0 


S (x)  * £ a u (x) 

nr  n nv  ’ 

n«r7 


1 


and  u (x)  represents  the  eigenfunctions  of  the  differential  operator 

( Ref  5:60-61).  Thus,  the  limit  of  the  expansion  of  eigenfunctions 

0 \ 

is  contained  in  a space,  H in  which  iStfily  the  functions  themselves 
are  required  to  be  square  integrable,  outside  the  solution  space  for 

the  problem,  in  which  the  functions  along\rith  their  first  and 

B \ 

second  derivatives  must  be  square  integrable  and  wl’>ich  must  satisfy 

\ 

the  boundary  conditions  on  the  problem,  from  which  the';  eigenfunct ions 

t 

u (x)  are  taken.  \ 

n 'r. 

In  a similar  manner,  the  trial  functions,  u(x)  , need  only  lie 
in  , since  Equation  (40)  can  be  integrated  by  parts  to  give 
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I(Q) 


1 /Tc— 

2 J [W 


■)  ♦ (U(X)}‘ 


2f(x) 


a(X)l 


dx 


(46) 


Thus,  the  functional  to  be  minimized  is  expressed  in  terms  of  the 
first  derivative.  Any  function,  u , will  be  admissible  in  the  mini- 
mizing process  provided  that  it  can  be  expressed  as  a limit  in  the  mean 

_ 2 

of  a sequence  of  functions,  Uj^  , which  lie  in  HR  (Ref  9:11).  Limit 

in  the  mean  is  the  equivalent  of  Equation  (44)  where  f(x)  is  the  limit 

in  the  mean  of  the  sequence  Sm(x)  . It  should  be  noted  that,  while 

2 

the  functions  u^  are  required  to  lie  in  Hg  , since  only  mean  square 
convergence  is  asked,  the  limit  u will  necessarily  satisfy  only  the 
Dirichlet  conditions  and  lie  in  the  space  H1  . Thus,  there  is  the 
advantage  that  the  functions,  u(x)  , can  be  constructed  from 
piecewise  linear  functions. 

The  value  I(u)  of  the  functional  I , then,  is  a limiting  case 
for  the  sequence  of  values  I(u^)  • The  function  u = u which  minimizes 

this  functional,  can  be  shown  to  satisfy  any  natural  boundary  conditions. 


where 


£ and  any 

v(x)  in 

Hg  , u(x)  > L-  o lies  in  H* 

and 

I(u)  < 

I(u) 

(47) 

Ku(x))  = 

I(u(x)  ♦ 

5 v(x)) 

Lr  -T 

X 

I(u)  ♦ 

2 E / ♦ uv  - fv  dx 

Q L.dx  dx  J 

(48) 

♦ 

L r 

52  l[ 

p2  ♦ 

1 
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Since  £ can  be  positive  or  negative,  the  middle  term  on  the  right 
hand  side  of  Equation  (48)  must  be  identically  equal  to  zero.  Thus, 

lo  [a?sf  * uv  - tv]dx“>  («> 

If  the  left  hand  side  is  integrated  by  parts,  the  result  is 

I [a7  5T  * uv  - f'"I  dx 

O L“ 

L L 

- (g-v)~|  + / ["(  - — ) V ♦ uv  - fv~]  dx  (50) 

*-o  o *-•  dx2 

The  right  hand  side  of  Equation  (50)  is  equal  to  zero  for  any  v(x)  in 

if  the  differential  equation  itself  is  satisfied  and  if,  assuming 

only  natural  conditions  were  originally  imposed  on  u(x)  , du/dx  = 0 

at  the  boundaries.  It  is  also  easy  to  see  why  any  essential  conditions 

must  be  imposed  on  v(x)  and  therefore  upon  u(x)  = u(x)  + £v(x)  , 

since  in  the  absence  of  a derivative  boundary  condition,  the  essential 

conditions  imposed  on  v will  be  needed  to  force  the  first  term  on 

the  right  of  Equation  (50)  to  vanish. 

N 

In  Figure  8,  0 (x)  , which  is  a piecewise  linear  function  defined 

in  terms  of  the  nodal  values,  can  be  substituted  into  I , and  I in 
turn  can  be  written  as  a sum  of  N functionals,  each  defined  over  a 
different  portion  of  the  domain.  These  portions  are  termed  "elements." 

That  portion  of  I defined  on  the  ith  element,  1^  , can  be 
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minimized  by  taking  the  derivative  of  1^  with  respect  to  each 
nodal  value  and  setting  the  result  equal  to  zero.  The  result  is  a 
system  of  algebraic  equations  which  must  be  solved  to  obtain  the  nodal 
temperatures.  Since  the  functional,  I , is  defined  over  the  entire 
interval,  the  minimization  process  takes  place  in  a mean  square  sense. 
This  is  an  important  distinction  between  finite-differences  and 
finite  elements.  In  the  finite-element  method,  the  errors  are 
distributed  over  the  entire  domain  and  the  solution  can  be,  but  is  not 
necessarily,  more  accurate  at  points  between  the  nodes  than  it  is  at 
the  nodes,  as  shown  in  Figure  10. 


11- 


Figure  10.  Example  of  a Finite-Element  Solution. 
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Application  of  the  Method.  Now  that  the  finite-element  method 
has  heen  introduced,  the  method  will  be  applied  to  the  heat  conduct* on 
equation. 

What  is  known  as  the  finite-element  method  is  actually  a group 
of  methods.  The  Rayleigh-Rltz  method  makes  use  of  the  associated 
variationul  principle  as  discussed  in  the  previous  section.  There  is 
another  method,  called  the  method  of  weighted  residuals  or  the  method 
of  Faedo-Galerkln,  which  is  nothing  more  than  the  discretization  in 
time  of  the  weak  form  of  the  problem: 

((§£  . iiH),  v)  - 0 (Si) 

36  3x2 

where  it  should  be  recalled  that  a and  x are  the  normalized  time 
and  space  variables.  The  weak  form  of  the  solution  does  not  require 
that  the  associated  functional,  I , he  minimized,  but  rather  that 
the  first  variation  of  I must  vanish  (Ref  9:9),  That  is, 

35  ■ 0 <*« 

The  expression  used  in  the  Galerkin  method  is  derived  in 
Appendix  E.  The  method  that  will  be  used  here  is  based  on  a 
variational  principle  that  was  suggested  by  Meyer  (Ref  3:399). 

Meyer  was  somewhat  obscure  in  the  justification  for  his  variational 
statement.  Indeed,  if  the  approach  he  seems  to  suggest  is  used,  the 
only  conclusion  which  may  be  drawn  is  that  no  variational  statement 


exists.  In  fact,  it  was  so  thought  by  Strang  and  Fix  (Ref  9:242)  and 
Washizu  (Ref  10) . It  turns  out  that  Gurtin  has  derived  a variational 
principle  for  the  heat  conduction  problem  which  involves  the  use  of 
convolutions  (Ref  11:255).  Meyer's  variational  principle  can  be 
justified,  however,  in  the  following  way.  If  the  space  and  time 
variables  are  separated,  then  the  heat  conduction  equation  can  be 
considered  to  be  the  equivalent  of  an  elliptic  equation  in  x at 
each  point  in  time,  and  a variational  statement  can  be  formulated  by 
neglecting  the  time  dependence.  The  variational  principle,  which  is 
simply  stated  by  Meyer,  is  derived  by  this  method  in  Appendix  E.  The 
resulting  functional,  which  is  to  be  minimized,  is 


n°)  ■ | /[^  * (4>2]  »* 


To  simplify  things,  the  integral  is  divided  into  two  integrals: 


1 f 8 ,, N 2 , 
7 J — (Q  ) dx 


o 36 


1 , 1 ,80N.2  , 

2 I dx 


The  following  formulation  is  developed  using  the  matrix 
notation  of  Meyer  (Ref  3:342-357). 
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I is  to  be  differentiated  with  respect  to  the  nodal 
temperatures  and  set  equal  to  zero  in  order  to  find  the  minimum 
within  the  space  SN  , If  uN  is  given  by 


then  the  minimizing  condition  is 


(56) 


(57) 


The  interval  [o,l]  is  divided  into  N elements,  as  shown  in 

Figure  11,  which  are  considered  individually.  Instead  of  differentiating 

N 

the  elemental  integrals  with  respect  to  each  component  of  u , a 
matrix,  , is  defined  by 


i 

i 


i 


Ji 
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and  used  as  follows: 


and 


dl, 


du 


N 


N m 
l D(l) 
i»l 


dl<« 


I £(1) 

i=l 


(59) 


(60) 


where  1^^  is  the  portion  of  I defined  on  the  ith  interval, 
(i,  i+1)  , and  where 


(61) 


The  temperature  distribution  within  each  element  must  be 
assumed.  Since  the  trial  functions,  (iN  , must  be  chosen  from  a 
subspace,  SN  , of  Hg  , they  must  be  continuous.  Piecewise  constant 
functions  cannot  be  used.  The  next  simplest  choice  is  piecewise 
linear  functions.  Therefore. 

0<4>  - cxW  ♦ c™  x (62) 

or 
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(63) 


T 

E £ 


(i) 


where 


£T  - [1  *] 


(64) 


and 


(65) 


The  coefficients  in  Equation  (62)  can  be  eliminated  by  considering  the 
set  of  equations  formed  by  substituting  Equation  (62)  into  Equation  (61) 
and  eliminating  c^*^  and  c^*^  • The  resulf  is 


u^ 


ET  i(ii  £(i) 


(66) 


where 


(67) 
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and 


6 xt  ■ XU1  - X. 


(68) 


The  subscript  on  Ax  can  be  dropped  since  the  interval  spacing  is 


assumed  to  be  constant.  Next,  the  derivative  of  u^  is  taken  with 


respect  to  x to  give 


k(i) 


3QV“'  = plR(i>u(i> 


dX  L*~  ~ 


(69) 


where 


(PT)  = t°  « 


(70) 


(i) 


Equation  (69)  can  then  be  substituted  into  1^  1 to  give 


i+1 


l{«  - \ / (PlR(i)  u(i))2dx 


(71) 


I.  is  then  differentiated  with  respect  to  uv  J to  give 


(i) 


dl 


(i) 


i+1  fnT  R(i)  U(i)l 

(px  ?.  y.  ] 


du 


d (Pi  R(l)  u(l))  dx 


(72) 


M 
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or 


Sot 


I1"  (p l r<‘>  a(l))  (px  r<»)t  ax 


Since  (p  R^  u^)  is  a scalar,  its  order  is  unimportant  and 


Equation  (73)  can  be  rearranged  as 


du^ 


J 1+1  (Px  R(l))T  (Px  R(l)  U(l))  dx 


And  since 


(A  B)T  > BT  AT 


Equation  (74)  is  equivalent  to 


Xi*‘  r<»T  p lp  T R<‘>  a'1')  dx 
*x  *x 


The  matrices  11^  and  are  independent  of  x and  can  be 

removed  from  the  integral: 
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JP* 


r'‘)T  [/  Xui  p,  £dx]  R(i)  0(i: 

— r — 


or,  if  the  integration  and  the  matrix  multiplication  are  carried  out 


11^  1 -1 
“l  1 1 .(i) 

TTI)  = Ax  \-\  lj  2. 


If  K''  ■*  , referred  to  in  the  literature  as  the  element  stiffness 


matrix,  is  defined  by 


KU)  a 


j_  p -r 
ax  L-i  i. 


then.  Equation  (59)  can  be  written 


dl,  N 


h - 1 o«  jc(i) 


“ N " ^ 

du  i=l 


Bseem 


L 


^ ■ l £(i)  K»>T  i 

dQ  i*l 


(81) 


where  was  defined  by  Equation  (56) . If  1C  , the  global  stiffness 

matrix,  is  defined  by 


N 


»(i)  „(i)  n(i)T 


K S l £'■ 

i*l 


(82) 


and  substituted  into  Equation  (81),  then 


dIl 


du 


K u 


(83) 


A similar  treatment  is  given  to  Equation  (60).  First,  since 


r(i)  _ 1 r *i*l  3 ra(i)>2 


iJ  w*caw)adx 


(84) 


X. 

1 


Leibnitz'  rule  for  differentiation  of  integrals  can  be  used  to  give 


«»>  - **/  *“.»“>)** 


(8S) 
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mm 


or 


1 

2 


_d_  f Xi*l 

de  1 


(PT  R(l)  Q(l))2  dx 


(86) 


The  derivative  of  this  expression  is  taken  with  respect  to 
to  yield 


di'1’ 

5F* 


d Ax 
d0  6 


2 

1 


(87) 


The  matrix 


M(i)  A Ax 

- T 


(88) 


is  referred  to  in  the  finite-element  literature  as  an  elemental  mass 
matrix.  This  definition  and  Equation  (87)  are  then  substituted  into 
Equation  (60) : 


dl 

du 


i 2(i)  4 (M'1’  ■ I o'1'  * (M<‘>  D<«T  QN)  (89) 

i-1  ao  i-i~  de  - - 
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Minimization  (actually,  the  second  derivative  would  have  to  be 
checked  to  prove  that  minimization  occurs)  is  forced  by  setting  the 
derivatives  in  Equations  (83)  and  (92)  equal  to  zero: 


t< 


< I 


The  tilde  has  been  dropped  since  the  solution  of  this  set  of  equations 

yields  the  closest  approximation  in  S^1  to  u . 

Equation  (94)  represents  a system  of  ordinary  differential 

equations.  These  equations  are  written  in  terms  of  the  nodal 

temperatures,  but  it  must  be  remembered  that  the  solution  is  assumed 

N 

to  vary  linearly  between  these  values,  for  it  is  u (x)  which  is 

N 

being  forced  close  to  u(x,0)  , not  just  the  nodal  values  of  u 

The  system  of  equations  (94)  could  be  solved  directly  for  a 
given  number  of  elements,  N . If  N is  large,  the  system  must  be 
solved  using  finite-differences  to  approximate  the  time  derivative, 
or  one  of  the  methods  suggested  by  Strang  and  Fix  may  be  used 
(Ref  9:244).  Three  common  methods  which  use  finite-differences  are 
the  Euler  method,  the  Crank-Nicolson  method,  and  the  fully  implicit 
method.  In  the  Euler  method,  forward  differences  are  used  to 
approximate  the  time  derivative  in  Equation  (94)  as 


, N k+1 
(u  ) 


A0 


(95) 


where  the  superscript,  k , indicates  the  present  time  level  and 
k>l  indicates  the  next  time  level,  after  a lapse  of  A0  • It  is 
important  to  note  that,  while  Equation  (95)  defines  (uN)lc^1  explicitly 
in  terms  of  temperatures  on  the  previous  time  level,  the  resulting 
scheme  after  substitution  of  Equation  (95)  into  Equation  (94)  is  not 
truly  explicit: 


1 

I 


i 
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M (uN)k+1  « (M  - K AO)  (uN)k 


(96) 


This  conclusion  is  a result  of  the  fact  that  M cannot  be  inverted 
without  destroying  its  tridiagonal  character.  Thus,  this  "explicit" 
scheme  has  no  advantage  over  implicit  schemes  in  the  number  of 
operations  required  to  obtain  a solution. 

A second  method  is  the  Crank-Nicolson  scheme,  which  approximates 
the  time  derivative  in  liquation  (94)  as  a central  difference  at  the 
point  k+1/2: 


d , N.k+4-  1,  d , N, 

3*r  (u  ) ^ - y (30^  ) 


N,  k + 1 


^N)k) 


, N k»l 

1m_1 


, N,  k 

;Ju 


AO 


(97) 


Substitution  of  Equation  (97)  into  Equation  (94)  gives 


1 


i 


(M  ♦ K 


AO, 


r N.k  + 1 

(U  ) 


,M  - K AO,  , N I 
(-  - y)  (U  ) 


(98) 


In  the  third  method,  a backward  difference  in  time,  at 
time  level  k*l  , gives 

i ( N k.l  . tuT"  - (99) 

3o  ' AO 

With  this  approximation,  Equation  (94)  becomes 
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(100) 


^ t N.k+1  ..  , N.k 

(M  ♦ K A9)  (u  ) = M (u  ) 

A general  scheme  which  can  generate  any  of  these  methods  is 
given  by 


(M  ♦ K aA0)  (uN)k>1  = (M  - K(l-a)A0)  (uN)k 

The  system  of  equations,  (101),  will  satisfy  any  natural 
boundary  conditions  imposed  upon  it,  automatically.  However,  all 
essential  conditions  must  be  forced.  In  the  case  of  problem  one, 
the  normalized  surface  temperature,  u^  , is  equal  to  zero  except 
for  the  initial  point  in  time,  at  which  a step  change  in  temperature 
occurs.  The  first  and  last  equation  in  the  system  (101)  can  be 
modified  as  follows:  The  matrices  and  l*  are  defined  by 

A * M + K aA0 


and 


B - M - K (l-a)A0 

These  matrices  are  substituted  into  Equation  (101),  then  the 
parameter  p , called  the  Fourier  modulus,  is  defined  by 


(101) 


(102) 


(103) 
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Because  of  the  constant  temperature  condition  on  the  boundaries 
the  matrices  become 


1 

0 

0 

0 

0 

U1 

0 

a22 

a23 

0 

0 

U2 

0 

a32 

a33 

a34 

0 

U3 

0 

0 

a43 

a44 

0 

U4 

0 

0 

0 

0 

1 

_U5  _ 

1 

0 

0 

0 

0 

UB 

0 

0 

0 

b22 

b23 

0 

0 

U2 

b21  UB 

a21  UB 

0 

b32 

b33 

b34 

0 

U3 

♦ 

0 

- 

0 

0 

0 

b43 

b44 

0 

U4 

b4S  UB 

a45  UB 

0 

0 

0 

0 

1 

_UB_ 

0 

0 

(106) 


where  the  lower  case  letters  represent  the  value  at  the  indicated 
position  in  A or  B and  ug  is  the  imposed  boundary  temperature. 

Since,  except  for  the  first  instant  in  time,  uR  = 0 , the  two 

column  matrices  on  the  right  hand  side  of  Equation  (106)  may  be  dropped. 

For  the  second  problem,  the  above  modification  is  needed  only 
for  the  left  boundary  condition,  since  the  Neumann  condition  on  the 
right  is  a natural  condition  for  the  method  and  is  satisfied  automatically. 
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Optimum  Implicit  Condition.  For  this  system  of  equations,  there 
is  an  optimum  choice  for  the  parameter,  a , shown  in  Figure  12.  If 
a is  chosen  according  to  the  formula  (derived  in  Appendix  D) 

* * T (1  * & 

then  the  resulting  expression  in  Equation  (106)  is  fourth  order 
accurate  at  the  nodes.  That  is,  the  truncation  error  at  the  nodes  is 
proportional  to  (Ax)4  . The  Euler,  Crank-Nicolson,  and  fully 

implicit  schemes  are  only  second  order  accurate. 

Error  Analysis 

General . One  of  the  objectives  of  this  project  was  to  compare 
the  accuracy  of  the  Crank-Nicolson  version  of  the  finite-element 
method  with  its  counterpart  in  finite-differences.  There  are, 
however,  several  fundamental  differences  between  finite-elements 
and  finite-differences.  First,  in  the  finite-element  method,  there 
are  two  ways  to  improve  the  accuracy  of  the  approximate  solution. 

The  first  is  to  decrease  the  size  of  the  interval  between  the  nodal 
points.  Ax  . This  procedure  has  a counterpart  in  finite-differences 
In  fact,  in  finite  differences,  if  the  limit  as  Ax  -*■  0 is  taken 
for  the  difference  equation,  and  p = A0/(Ax)2  , then  the  difference 

equation  should  converge  to  the  differential  equation  in  the  limit. 

For  the  finite-element  method,  there  is  an  additional  procedure 
for  improving  the  accuracy  of  the  solution.  Since  the  temperature 
distribution  is  assumed  to  have  a certain  shape  within  each  element, 
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in  the  case  of  a polynomial  approximation,  a higher  order  polynomial 
will  give  better  accuracy.  This  procedure  has  no  counterpart  in 
finite-differences,  in  which  only  the  nodes  are  considered. 

Secondly,  the  error  in  the  finite-element  method  is  most 
naturally  measured  over  the  entire  interval  rather  than  in  the 
pointwise  sense.  The  method  was  intended  to  minimize  the  error 
between  the  approximate  solution  and  the  exact  solution  in  a mean 
square  sense.  This  fact  must  be  considered  when  comparing  the  error 
in  the  two  methods. 

In  order  to  compare  linear  finite-elements  to  ordinary  finite- 
differences,  a common  yardstick  must  be  used.  Three  error  norms 
were  chosen  which  could  be  applied  to  both  methods,  and  one  is 
also  used  which  applies  only  to  finite-elements. 

Finite-Difference  Error  Analysis.  The  first  measure  of  the 
error  for  finite-differences  is  the  pointwise  error.  The  estimate 
for  this  error  is  given  by  a truncation  error  analysis.  Actually, 
the  pointwise  error  is  a sum  of  two  types  of  error.  The  total 
pointwise  error  will  be  referred  to  as  the  discretization  error. 

This  error  is  composed  of  round-off  error,  which  results  from  carrying 
only  a finite  number  of  significant  figures  in  all  of  the  computations, 
and  truncation  error,  which  results  from  elimination  of  the  higher 
order  derivatives  when  Taylor's  series  is  used  to  approximate  a 
differential  equation.  Crandall  suggests  that  discretization  error 
can  be  decreased  by  using  a smaller  nodal  spacing,  Ax  , in  the 
difference  equation,  but  warns  that  if  the  time  interval  is  also 
decreased,  to  maintain  the  same  stability,  then  more  computations 
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will  be  required  to  cover  the  time  domain.  Thus,  when  the  number 
of  time  steps  is  large,  round-off  error,  while  usually  much  smaller 
in  magnitude  than  the  truncation  error,  could  become  significant. 

One  solution  to  this  problem  is  to  increase  the  number  of  significant 
figures  carried  in  the  computation  (Ref  12:170).  Thus,  while  the 
discretization  error  is  actually  what  is  being  measured,  it  is 
reasonable  to  consider  this  error  as  roughly  equivalent  to  the 
truncation  error.  The  truncation  error  is  derived  in  Appendix  C 
for  the  general  expression.  Equation  (32),  and  results  in  the 
following  estimate: 


et(iAx,kA0)  = Cj  [j  - ap  - -jj]  Ax2 
♦ 0 (ax)**  + ... 


(108) 


where 


c 


1 


32u 

302 

x * iAx 
0 * kA0 


and  0(Ax)4  indicates  that  the  first  term  which  is  neglected 
is  proportional  to  (Ax)4  . This  error  estimate  applies  to  all 

nodal  points  in  a Dirich let  problem.  In  the  secondary  problem, 
the  Neumann  condition 
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3u 

3x 


= 0 


x = 1 


(109) 


can  be  analyzed  by  using  the  type  of  truncation  error  analysis 
suggested  in  Appendix  C.  The  truncation  error  at  the  right  boundary 
for  Equation  (37)  is  thus  given  by 

e (NAx.kAG)  = (2ap  + 4)  Ax  — - 
* 3x3 

where  the  condition  3u/3x  =0  at  x = 1 must  be  applied  to  cancel 
the  lower  order  terms.  This  result  agrees  with  Crandall  (Ref  12:266) 
who  indicates  an  order  of  Ax  is  lost  for  each  degree  in  the 
derivative  boundary  condition. 

The  error  estimates  given  above  apply  to  all  of  the  error 
measures  used  for  the  finite-difference  method.  The  first  error 
measurement  is  called  pointwise  error.  This  is  simply  the  total 
discretization  error  at  a given  nodal  point.  It  is  defined  as  the 
difference  between  the  exact  analytical  solution  and  the  discrete 
approximate  solution  at  that  point.  If  the  round-off  error  is  assumed 
to  be  negligible,  this  error  is  equivalent  to  the  truncation  error. 

The  second  error  measure  is  defined  to  be  the  maximum  error 
between  the  exact  solution  and  the  finite-difference  solution  taken 
at  any  node.  In  the  finite-element  literature,  the  continuous  analog 
of  this  error  measure  is  called  the  Tchebycheff  norm.  Because  it 


NAx 

1/AA 


(110) 


56 


is  defined  here  at  nodal  points  only,  it  is  estimated  by  the 
truncation  error  estimates  previously  derived. 

The  last  error  measure  for  finite-differences  is  called  here 
a generalized  mean  error.  Actually,  it  is  nothing  more  than  the  sum 
of  the  absolute  values  of  the  discretization  errors  at  each  of  the 
non-zero  nodes.  This  generalized  mean  was  devised  to  compare  the 
error  at  the  first  interior  node,  x = .1  , which  is  measured  in 
the  pointwise  sense,  and  the  error  at  all  the  nodes.  The  idea  here 
was  to  investigate  whether  or  not  a higher  or  lower  order  of 
convergence  would  be  seen  for  the  first  interior  node  compared  to 
the  convergence  at  all  of  the  nodes.  There  was  still  another  reason 
for  this  error  measure  which  will  be  discussed  in  the  next  section. 

In  Equation  (108),  it  is  interesting  to  note  the  order  of 
accuracy  for  each  of  the  four  finite-difference  schemes  which  have  been 
discussed:  explicit,  Crandall,  Crank-Nicolson,  and  fully-implicit. 

For  the  explicit  scheme,  a = 0 ; therefore 

e«.(iAx,kA0)  = & - ~)  (Ax)2—  + 0(Ax)4  +. . . (Ill) 

1 2 12  3x2 

x = iAx 
e = kA0 


or  second  order  accurate.  For  the  Crandall  method,  o = (l-l/(6p))/2  , 

and 
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And  last,  for  the  pure  implicit  method,  a = 1 , and  therefore, 

e (iAx,kA0)  = " T5O  (Ax)2  ~ ♦ 0(Ax)4  + . . • (114) 

1 u ax2 

x * iAx 
0 * kA0 

fe 

It  is  interesting  to  note  that,  with  the  exception  of  the  Crandall 
method  which  is  fourth  order  accurate,  all  of  the  above  schemes  are 
second  order  accurate  with  respect  to  truncation  error. 

The  next  section  will  deal  with  error  estimates  for 
finite-elements . 

Finite-Element  Error  Analysis.  As  in  the  finite-difference  method, 
the  first  error  measure  which  will  be  considered  will  be  the  pointwise 
error.  This  error  measure  is  alien  to  the  finite-element  method  because 
the  concept  behind  the  method  is  to  minimize  the  error  everywhere  in 
the  region,  not  just  at  the  nodes.  However,  error  in  the  pointwise 
sense  can  be  estimated  by  treating  the  individual  equations  in  the 
1 system  (106)  as  if  they  were  simple  difference  equations.  In 
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Appendix  D,  the  truncation  error  is  derived  for  the  ith  equation  in 
this  system.  The  result  is 


«t(iAx,kA0) 


l \ - <>p  + -nrl  (Ax)2  — 
1 12  3x2 


♦ 0(Ax)4 


x = iAx 
6 * kA0 


This  error  estimate  applies  to  all  nodes  for  a Dirichlet  problem. 

For  the  secondary  problem,  the  truncation  error  for  the 
Neumann  boundary  condition  is  given  by 


et(iAx,NA0) 


6op  Ax 


x = NAx 
0 = kA0 


(115) 


(116) 


This  order  of  accuracy  is  the  same  as  the  predicted  order  of  accuracy 
in  the  norm  for  a derivative  boundary  condition  (Ref  9:118). 

The  error  estimates  given  in  Equations  (115)  and  (116)  apply  to 
the  pointwise  error,  the  discrete  Tchebycheff  norm  (maximum  error  at 
any  node) , and  the  generalized  mean  error  (as  defined  for  the  finite- 
difference  method).  However,  as  mentioned  previously,  there  are  other 
methods  for  measuring  the  error  in  the  finite-element  method.  The 
method  used  here  is  called  the  L norm  or  sometimes  the  H°  norm. 

It  is  defined  as  follows  (Ref  9:5): 
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(117) 


I |u-uN| l0  = [/0  (U  - UN)2  dx  ] /2 

Since  the  exact  analytical  solution  for  u(x,0)  is  given  in  terms  of 
an  infinite  series,  some  special  questions  on  the  existence  of 
| |u-u^| |q  must  be  answered.  The  computation  is  involved;  therefore, 
it  has  been  placed  in  Appendix  G.  Strang  and  Fix  develop  a theorem 
which  bounds  the  error  in  the  norm  for  a piecewise  linear  finite- 

element  space  by  (Ref  9:250) 

||u-uN||o  < C (Ax)2  (118) 

where  the  constant  C is  constant  only  with  respect  to  the  spatial 
domain,  but  can  be  a different  value  at  each  point  in  time.  This 
order  of  accuracy  applies  only  when  the  full  admissible  space  of  functions, 
u , is  Hg  , and  , the  space  of  all  functions  uN  , is  a 

piecewise  linear  subspace  of  Hg  • 

Stability  Analysis 

General.  Since  only  a finite  number  of  significant  figures  can 
be  carried  out  by  the  computer  in  a calculation,  every  time  a 
calculation  is  performed,  there  is  the  chance  of  introducing  an  error. 

This  error  is  called  round-off  error.  If  a series  of  finite-difference 
or  finite-element  computations  was  carried  out  using  an  infinite  number 
of  significant  figures,  the  result  would  be  the  exact  solution  of  the 
set  of  equations.  If  the  magnitude  of  the  difference  between  this 
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exact  numerical  solution  and  the  truncated  numerical  solution,  which 
would  be  generated  by  a computer  using  a finite  number  of  significant 
figures,  grows  exponentially  as  the  calculation  proceeds,  then  the 
numerical  scheme  is  termed  unstable. 

There  are  several  methods  for  treating  stability,  but  by  far 
the  best  method  is  one  which  deals  with  the  complete  numerical  scheme 
including  the  boundary  conditions.  The  method  used  here  is  an 
adaptation  of  one  used  by  Crandall  (Ref  12:382).  Crandall,  however, 
indirectly  addresses  this  question  using  the  method  of  separation  of 
variables  to  isolate  the  spatial  dependence  from  the  temporal  in  the 
case  of  a parabolic  problem.  He  then  analyzes  the  eigenvalues  of  the 
separated  spatial  system  of  equations.  An  eigenvalue  greater  than  one 
in  value  will  cause  steady,  unbounded  growth  in  its  associated  spatial 
eigenfunction,  assuming  that  eigenfunction  was  excited  by  the  initial 
conditions  of  the  problem  or  was  introduced  by  round-off  errors  during 
the  computation.  In  other  words,  the  spatial  mode  associated  with  the 
eigenvalue  will  be  amplified  and  be  of  the  same  sign  after  each 
computation.  If  the  value  of  an  eigenvalue  lies  between  zero  and  one, 
it  will  cause  steady  decay  of  the  corresponding  spatial  mode.  The 
amplitude  of  that  mode  will  be  smaller  in  magnitude  and  of  the  same 
sign  after  each  computation.  If  an  eigenvalue  has  a magnitude  between 
zero  and  minus  one,  its  associated  eigenfunction  will  be  smaller  in 
amplitude,  but  alternate  in  sign  with  each  step  in  the  computation. 
Under  these  conditions,  the  solution  is  said  to  undergo  stable 
oscillations.  Finally,  if  an  eigenvalue  has  a value  less  than  minus 
one,  then  its  eigenfunction  will  undergo  unstable  oscillations  where 
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the  amplitude  will  increase  in  magnitude,  but  have  an  alternating 
sign  with  each  step. 

Here,  the  error  in  the  solution  is  the  thing  of  interest. 

Round-off  errors  should  not  grow  exponentially  for  a stable  solution. 

If  Equations  (32)  and  (106)  are  written  in  terms  of  error  vectors  where 
the  vector  represents  an  error  introduced  by,  say,  round-off, 

and  e^  represents  the  new  error  vector  after  solution  of  the  set  of 
equations,  then  the  following  equation  can  be  written: 


A 


B e 


(119) 


or 


e.  » C e 
—1  c 


(120) 


where  £ * (A  1 B)  . If  e^  is  expanded  in  terms  of  the  eigenvectors, 
, of  C then 


n 

£ I 'i , 

i»l  *i 


where  c^  is  a constant  and  4k  is  the  ith  eigenvector  of  £ 
the  definition  of  an  eigenvalue,  Equation  (119)  can  be  written 


•i  ■ ci  ‘i 
— i-1  — 


(121) 


With 


(122) 
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where  X^  is  the  ith  eigenvalue  of  the  matrix  C . Similarly 


®2  * — el 


l ci  Xi 
i-1  — 


(123) 


and  after  k computations 


ek  ■ I; ci  *i 

— ia  1 — 


(124) 


This  demonstrates  that  the  values  of  the  eignevalues  of  the  iteration 
matrix,  C , determine  the  growth  or  decay  of  errors  just  as  the 
eigenvalues  of  the  separated  spatial  system  did  before.  Separating 
the  variables  in  the  difference  equation  is  not  necessary  in  this  case, 
however.  All  that  is  required  is  to  solve  for  the  eigenvalues  of  the 
iteration  matrix  £ . 

The  stability  analysis,  as  discussed  here,  will  be  applied 
only  to  the  primary  problem. 


Analysis  of  the  Finite-Difference  Formulation.  Equation  (32)  can 


be  written  in  vector  notation  as 


=>  B u 


(12S) 


Both  A and  B are  tridiagonal  of  the  form 


( 


b a 

aba 

0 0 0 

aba 
a b 


(126) 


Smith  (Ref  8:65)  gives  the  eigenvalues  of  a matrix  of  this  type  as 


c 


A = b + 2 a cos  C^V)  , n = 1,2,...,N 
n N+l' 


where  N is  the  order  of  the  matrix.  The  iteration  matrix, 
C * A~*  B , in  Equation  (125)  has  eigenvalues 


(A  ) 

C n 


<Vn 

<Vn 


or  from  Equation  (32) 


<Vn 


(1+2  pa)  + 2 (-pa)  cos  nn 

N+l 


(1  ♦ 2p(l-a))  + 2(p(l-a))  cos  rnr 

N+l 


(127) 


(128) 


(129) 


The  analysis  will  be  done  for  the  limiting  case  where  N ■+■  •.  , that 
is,  for  a large  system  of  equations.  The  eigenvalues  for  a system 
of  40  equations  are  very  close  to  the  limiting  values.  The  eigenvalues 
of  such  a system  are  bounded  above  by  one.  The  thing  of  interest  here, 
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then,  is  the  minimum  eigenvalue,  which  is  given,  for  a given  p and 
a given  a , by  Equation  (129)  when  n = N . In  the  limit,  then 


lim  , Nir 

N-~  C0SW 


-1 


(130) 


Thus,  the  eigenvalue  of  interest  is  given  by 


cv- 


1 + 4pq 
1 - 4p(l-a) 


(131) 


which  will  be  termed  the  critical  eigenvalue  and  where  the  subscript 
indicates  that  n = N = » . Since  a gives  the  "degree  of  implicitness" 

of  Equation  (32),  the  entire  family  of  finite-difference  expressions 
can  be  investigated.  Figure  13  shows  graphically  the  results  of 
this  analysis  for  the  pure-implicit,  Crank-Nicolson,  Crandall,  and 
explicit  formulations.  Table  1 gives  the  values  for  the  Fourier 
modulus,  p , which  will  cause  oscillation  or  instability  in  each 
scheme.  The  values  given  here  are  the  same  as  those  given  by 
Crandall  using  the  separation  of  variables  technique  (Ref  2:319). 
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CRITICRL  EIGENVALUE 

rl.OO  -0.71  -0.43  -0.14  0.14  0.43  0.71 


Figure  13.  Stability  Curves  for  Finite-Differences. 
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TABLE  I 


Oscillation  and  Instability  Limits  for 
the  Fourier  Modulus  in 
the  Finite-Difference  Formulation. 


Limits 

Explicit 

Crandall 

Crank - 
Nicolson 

Pure- 

Implicit 

Oscillation  Limit, 

No 

p < X 

for  no  oscillation 

0 

.25 

.3333 

.5 

Oscillations 

Stability  Limit, 

p < X 

for  stable  scheme 

1 

.5 

Always 

Stable 

Always 

Stable 

Always 

Stable 

In  this  formulation,  the  first  and  last  equations  in  this  system  can 
be  dropped.  The  number  of  unknowns  can  be  reduced  by  two  by 
eliminating  the  known  boundary  values,  uD  . All  of  these  changes  can 

D 

be  made  without  altering  the  numerical  solution.  If  Equation  (132) 
represents  the  new  system  after  these  changes  have  been  made,  then  the 
resulting  system  (132)  can  be  analyzed  exactly  as  was  done  for  the 
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finite-difference  case.  Again,  both  A and  11  are  tridiagonal  of  the 
same  form  as  Equation  (126) , If  the  same  rationale  which  was  previously 
used  in  the  finite-difference  case  is  applied  here,  then  the  critical 
eigenvalue  of  the  iteration  matrix  = A-*  is  given  by 


(V» 


1 + 12  pa 
1 - 12p(l-a) 


(133) 


Figure  14  shows  the  relationship  between  the  various  finite-element 
schemes.  Table  II  gives  the  oscillation  and  instability  limits  in  terms 
of  the  Fourier  modulus,  p . A discussion  of  the  stability  curves  and 
limits  will  be  reserved  for  the  results  section. 

TABLE  II 

Oscillation  and  Instability  Limits  for 
the  Fourier  Modulus  in  the 
Finite-Element  Method. 


Limits 

Euler 

Crank- 

Nicolson 

Optimum 

Implicit 

Pure 

Implicit 

Oscillation 

Limit,  p < x 

for  no 

oscillation 

.08333 

. 16667 

.33333 

Never 

Oscillates 

Stability 

Limit,  p < x 

for  a 

stable 

scheme 

. 16667 

Always 

Stable 

Always 

Stable 

Always 

Stable 
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CRITICAL  EIGENVALUE 

,-i .00  -0.71  -0.43  -0.  U 0.14  0.43  0.7i 


Figure  14.  Stability  Curves  for  Finite  Elements  . 
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I 

III , Procedure 

Approach 

Initial  Phase.  The  first  step  in  this  investigation  began  with 
the  programming  of  the  finite-difference  method.  This  phase  was 
relatively  easy  since  the  necessary  background  for  this  method  had  been 
obtained  in  previous  course  work.  Concurrently,  a literature  search 
was  initiated  along  with  a study  of  conduction  heat  transfer,  and  a 
study  of  the  basic  principles  of  the  finite-element  method.  Meyers 
(Ref  3)  was  the  text  for  the  initial  study  of  the  finite-element 
method.  His  treatment  was  on  an  introductory  level;  a more  complete 
treatment  of  the  calculus  of  variations  was  found  in  Forray  (Ref  14) . 

Second  Phase.  This  phase  began  with  the  programming  of  the 
finite-element  method.  The  study  of  the  finite-element  method  was 
continued  during  this  phase,  as  indeed  it  continued  throughout  the 
entire  project.  Initial  results  collected  during  this  phase  were 
encouraging  in  that  they  indicated  that  there  may  be  a point  of  minimum 
error  for  the  method.  The  difficulty  lay  in  how  to  predict  this 
condition.  The  solution  of  the  problem  is  a function  of  x , 

0 , p , Ax  , A0  , and  a . Not  only  that,  but  the  error 
could  be  measured  in  a number  of  ways.  This  meant  that  there  was 
an  enormous  amount  of  data  to  be  analyzed  and  a number  of  ways  to 
analyze  it.  As  the  sifting  of  this  data  took  place,  various  methods 
of  presentation  of  the  data  were  considered.  It  was  decided  that  a 
method  used  by  Campbell,  Kaplan,  and  Moore  (Ref  4:325)  would  provide 

( 
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the  most  information  concerning  convergence  of  the  numerical  solution 
to  the  true  solution  as  Ax  becomes  smaller.  It  should  be  mentioned 
that,  as  Ax  gets  very  small,  round-off  error  begins  to  dominate 
the  solution  such  that  there  exists  a value  of  Ax  , for  a given 
number  of  significant  figures  carried  in  calculation,  beyond  which 
the  approximate  solution  gets  less  accurate.  Campbell,  et.  al., 
defined  the  discretization  error  ratio  as  the  ratio  of  the  error 
incurred  when  one  subdivision  of  the  space  domain  is  used  to  the  error 
at  the  same  point  when  the  number  of  nodes  has  been  doubled.  For  the 
1*2  norm  in  finite-elements,  it  is  somewhat  of  a misnomer  since  the 
assumed  solution  is  not  discrete,  but  continuous;  however,  the 
terminology  is  retained  and  is  understood  to  be  the  error  ratio  over 
the  entire  spatial  domain  after  the  number  of  elements  has  been 
doubled.  At  about  this  point  in  time,  it  was  discovered  that  by 
making  a simple  truncation  error  analysis  of  the  finite-element  equations, 
as  if  they  were  simple  difference  equations,  and  setting  the  coefficient 
of  the  0(Ax)2  term  in  the  resulting  expression  equal  to  zero,  that 
an  0(Ax)4  scheme  could  be  obtained.  This  was  very  significant  since 

the  optimum  value  of  the  parameter  a could  be  predicted  rather  than 

estimated  from  empirical  results.  Indeed,  such  a large  amount  of 
data  would  be  needed  to  find  the  optimum  condition  as  a function  of 

all  the  variables  in  the  problem,  that  the  search  was  sure  to  be  very 

time  consuming  and  perhaps  even  inconclusive.  It  turns  out  that  the 
optimum  value  of  a is  a function  of  the  Fourier  modulus,  p , only, 
and  is  invariant  with  respect  to  the  other  five  variables,  at  least 
in  theory. 
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Third  Phase.  The  third  phase  of  the  investigation  began  with 
a study  of  error  in  both  the  finite-difference  and  finite-element 
methods.  More  emphasis  was  placed  on  finite-elements  for  which 
the  error  norms  are  more  involved.  Next,  the  question  of  stability 
was  investigated.  Of  the  methods  considered.  Von  Neumann,  Dusinberre, 
and  matrix  eigenvalue,  one  method  stood  out  as  being  straightforward 
and,  in  fact,  more  precise.  This  was  the  matrix  eigenvalue  method. 

Of  course,  the  other  methods  are  useful  if  the  eigenvalues  of  the 
iteration  matrix  are  too  difficult  to  obtain. 

Also  during  this  phase,  the  results  of  the  primary  problem  were 
analyzed.  Once  the  optimum  value  of  a was  determined  by  theory, 
computer  calculations,  which  investigated  the  behavior  of  the  solution 
in  the  neighborhood  of  the  optimum  value  of  a , were  begun.  The 
initial  results  of  this  investigation  were  disappointing  since  they 
showed  little  improvement  in  the  discretization  error  ratio,  or 
rate  of  convergence,  for  the  predicted  optimum  value  of  a . The 
source  of  the  difficulty,  it  was  discovered,  was  the  discontinuity 
between  the  initial  condition  and  the  boundary  conditions.  This 
problem  is  discussed  by  Smith  (Ref  8:48-49)  and  by  Campbell,  Kaplan, 
and  Moore  (Ref  4:325-326).  Smith  suggests  two  methods  for  handling 
this  problem.  The  first  involves  a transformation  of  the  independent 
variables  from  (x,0)  to  (X,0)  where 


X * x(0) 


(134) 


. .*/» 


0-0 


(135) 
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The  result  of  this  transformation  is  an  expansion  of  the  origin 
(0,0)  onto  the  positive  side  of  the  new  X axis  while  the  old  x 
axis  is  concentrated  at  a point  at  infinity  on  the  X axis.  The 
jump  condition  at  the  origin  has  been  transformed  into  a continuous 
change  defined  over  the  positive  side  of  the  X axis.  Smith  suggests, 
also,  that  an  alternative  would  be  to  calculate  an  analytical  solution 
which  is  continuous  in  the  neighborhood  of  the  jump  condition.  A 
third  approach,  subdividing  the  space  time  grid  for  the  first  time 
step  as  suggested  by  Campbell,  et.  al.,  is  the  one  chosen  for  this 
investigation.  This  last  method  was  chosen  because  of  the  ease 
with  which  it  could  be  incorporated  into  the  computer  programs  written 
during  phases  one  and  two. 

It  was  felt  that  the  merits  of  this  investigation  lay  not  in 
overcoming  the  problem  with  the  discontinuity,  but  rather  in  predicting 
and  verifying  the  existence  of  an  optimum  value  of  the  parameter  o 
for  the  finite-element  method.  For  this  reason,  much  of  the  investi- 
gation was  concerned  with  a modification  of  both  the  primary  and 
secondary  problems.  The  modification  was  made  by  simply  substituting 
the  exact  analytical  solution  after  the  first  time  step.  In  effect, 
this  transformed  the  original  problem  into  a new  problem  in  which 
no  discontinuity  existed  between  the  initial  condition  and  the  boundary 
conditions.  After  this  transformation,  the  discretization  error  ratio 
was  found  to  approach  the  predicted  values  at  the  nodes. 

There  was  another  method  which  was  used  in  an  attempt  to  reduce 
the  effect  of  the  discontinuity.  It  was  pointed  out  that,  since  the 
constant  initial  temperature  in  problem  one  was  represented  by  setting 
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Ug  equal  to  zero  in  Equation  (106),  that  this  method  underestimates 
the  initial  temperature  distribution.  Since  the  coefficients  of 
b2j  and  b^,.  in  the  first  column  matrix  in  Equation  (106)  should 
more  properly  be  equal  to  unity  for  the  first  time  step,  it  was  thought 
that  a calculation  made  using  that  modification  would  give  better 
convergence  rates  that  the  original  version.  This  modification  was 
used  only  for  the  first  time  step  after  which  the  original  scheme 
was  used.  The  results  are  shown  in  Figures  15  and  16  for  two 
different  points  in  time.  It  can  be  seen  that  the  effect  is  to 
overestimate  the  temperature  distribution.  This  modification  was  not 
pursued  any  further. 

Fourth  Phase.  The  last  phase  of  the  project  began  with  a 
reprogramming  of  the  two  methods  to  solve  the  secondary  problem.  After 
that  was  completed,  the  next  step,  of  course,  was  to  analyze  the  results 
from  the  second  problem  for  both  finite-elements  and  finite-differences. 
The  same  types  of  problems  which  were  encountered  with  the  primary 
problem  were  also  encountered  with  the  secondary  problem  and  were 
handled  in  the  same  way. 

Computer  System  and  Programs 

Computer.  The  computer  system  used  for  this  project  is  one 
designed  by  the  Control  Data  Corporation.  It  consists  of  input  and 
output  devices,  peripheral  processors,  and  two  central  processors. 

The  central  processors  are  a CDC  6613  and  a CDC  CYBER  74  which 
operate  in  parallel.  Each  has  131,000  60-bit  words  of  central  memory. 
Magnetic  disc  and  drum  storage  were  used  as  temporary  storage  devices. 
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SPLINE  FIT  OF  EX.  80LN.  41  PTS 


MODIFIED  I.C.  VS.  ORIGINAL 
ALPHA=. 5.  TIME=. 01 


-Original  Method 

-Exact  Solution 

-Method  Using  the 
Modified  Initial 
Condition 


IS  rh«  Exact  Analytical  Solution  Compared  to  the  Numerical  Solution 
'•in*  'nginal  and  Modified  Approximations  for  the  Initial  Conditions 
•*  • >i  for  the  Crank -Nicolson  Version  of  the  Finite-Element  Formulation. 
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MODIFIED  I.C.  VS.  ORIGINAL 
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-Original  Method's 
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-Method  Using  the 
Modified  Initial 
Condition 


^Too  0.20  0.40  0.60 

X (NORMALIZED) 


0.60 


1.00 


Fig.  16.  The  Exact  Analytical  Solution  Compared  to  the  Numerical  Solution 
Using  the  Original  and  Modified  Approximations  for  the  Initial  Conditions 
at  0 ■ .02  for  the  Crank-Nicolson  Version  of  the  Finite-Element  Formulation. 


Computer  Programs.  Four  major  programs  were  written  during  this 
investigation.  TVo  finite-difference  programs  and  two  finite-element 
programs . There  were  also  numerous  other  programs  which  were  written 
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for  stability  analysis  and  plotting  purposes.  The  programs  for  the 
secondary  problem  were  primarily  modifications  of  the  ones  written 
for  the  primary  problem.  The  language  which  was  used  for  all 
the  programs  was  FORTRAN  EXTENDED. 

Cost  and  computer  run  time  are  certainly  of  interest  in  most 
computer  work,  but  because  the  thrust  of  this  investigation  was  in 
the  direction  of  investigating  the  behavior  and  accuracy  of  the 
finite-element  solution  as  a function  of  the  parameter  a , there 
was  no  effort  to  compare  cost  and  run  time  for  various  options.  The 
main  reason  for  this  was  due  to  the  fact  that  the  programs  which  were 
written  were  designed  for  generality  and  make  use  of  mass  data 
storage  on  temporary  files.  For  these  reasons,  an  accurate  comparison 
of  run  time  and  cost  could  not  be  made.  However,  except  for  the  explicit 
finite-difference  scheme,  all  of  the  methods  considered  were  implicit 
with  tridiagonal  matrixes  and  thus  should  have  been  roughly  comparable 
with  respect  to  computation  time.  Also,  modifications  would  need  to 
be  made  for  production  codes.  Certainly,  the  user  would  need  to 
store  only  two  columns  of  data  to  represent  the  symmetric  tridiagonal 
matrixes  which  are  generated  by  the  two  methods  as  used  here.  Flow 
charts  are  given  in  Figures  17  and  18  for  the  primary  problem 
using  finite-differences  and  finite-elements,  respectively. 
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Figure  18.  Finite-Element  Computer  Program  Flow  Chart. 
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IV.  Results 

Stability  Analysis 

The  results  of  the  stability  analysis  have  already  been  presented 
in  graphical  and  tabular  form  in  Chapter  II.  It  will  be  pointed  out 
here,  however,  that  it  is  very  interesting  to  note  that  the  stability 
curve  for  the  optimum  implicit  finite-difference  scheme  proposed  by 
Crandall  coincides  exactly  with  the  optimum  implicit  finite-element 
scheme  proposed  here.  Further,  it  is  interesting  to  note  that,  while 
the  Crandall  method  is  less  stable  than  the  Crank-Nicolson  method  for 
finite-differences,  the  optimum  implicit  method  suggested  for  finite- 
elements  is  more  stable  than  the  Crank-Nicolson  method  for  finite- 
elements.  In  general,  the  finite-element  method  is  more  restrictive 
with  respect  to  stability;  that  is,  the  oscillation  limits  and 
instability  limits  occur  for  smaller  values  of  the  Fourier  modulus, 
p - A0/(Ax)2  . 

Error  Analysis 

The  results  for  the  error  analysis  occupy  the  bulk  of  the 
discussion  in  this  chapter.  A very  large  number  of  plots  have  also 
been  generated  during  the  investigation  and  have  been  placed  in 
Appendix  H for  reference.  These  plots  are  offered  as  proof  for  the 
theories  set  forth  in  this  paper  and  only  the  trends  will  be  identified. 
Specific  comments  on  each  plot  will  not  necessarily  be  made. 

First,  one  of  the  objectives  of  this  research  was  to  compare 'the 
finite-element  method  with  the  finite-difference  method  with  respect 
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to  accuracy.  Tables  III  through  VI  give  the  results  of  this  analysis 
for  two  points  in  time  for  the  case  where  the  space  domain  is 
divided  into  20  intervals.  Surely,  a more  complete  picture  could  be 
presented,  and  the  data  is  available;  but  the  real  thrust  of  this 
project  was  in  verifying  the  new  theories  developed. 


Table  III 

Error  Comparisons  for  the  Various  Methods 
for  0 = .04  and  p = 1.0  in  the  Primary  Problem. 


Method* 

Pointwise 

Error,  x = . 1 

Maximum  Error 
at  Any  Node 

Generalized 
Mean  Error 

L2  Error 

Norm  (FE) 

CNM-FD 

CNM-FE 

7.7494  x 10~6 
1.7696  x 10-2 

1.4877  x 10-3 
2.6208  x 10*3 

6.0711  x 10-3 
1.8231  x 10"  2 

3.1957  x 10-3 

OIM-FD 

OIM-FE 

8.7106  x 10"4 
8.7106  x 10"4 

1.5579  x 10-3 
1.5579  x 10’3 

1.2018  x 10“2 
1.2018  x 10-2 

2.6011  x 10"3 

FIM-FD 

FIM-FE 

-5.60758  x 10‘3 
-3.6498  x 10"3 

7.1078  x 10~3 
4.5849  x 10-3 

3.6781  x lO-2 
2.4720  x lO"2 

2.5395  x lO'3 

*CJW  = Crank-Nicolson  Method 
OIM  * Optim  jn  Implicit  Method 
FIM  = Fully  Implicit  Method 
FD  * Finite-Differences 
FE  * Finite -Elements 


Table  IV 


Error  Comparisons  for  the  Various  Methods 
for  6 = .08  and  p = 1.0  in  the  Primary  Problem. 


Method* 

Pointwise 
Error,  x - .1 

Maximum  Error 
at  Any  Node 

Generalized 
Mean  Error 

L2  Error 

Norm  (FE) 

CNM-FD 

CNM-FE 

6.0493  x 10~5 
7.1111  x 10~3 

3.0606  x 10'4 
2.0911  x 10"3 

1.6801  x 10-3 
1.3683  x 10‘2 

2.3863  x 10-3 

OIM-FD 

OIM-FE 

3.8710  x 10"4 
3.8710  x 10‘4 

1.1952  x 10-3 
1.1952  x 10-3 

7.6779  x 10"3 
7.6779  x 10-3 

1.7356  x 10-3 

FIM-FD 

FIM-FE 

-1.9654  x 10“3 
-1.2764  x IQ"3 

4.8803  x 10-3 
3.1811  x IQ'3 

3.4197  x 10-2 
2.2256  x IQ-2 

1.6441  x 10-3 

♦See  Table  III 


Table  V 

Error  Comparisons  for  the  Various  Methods  for  0 = .04  and 
p * 1.0  in  the  Primary  Problem  and  where  the  Exact  Analytic  Solution 
has  been  substituted  at  the  First  Time  Step. 


Method* 

Pointwise 
Error,  x = .1 

Maximum  Error 
at  Any  Node 

Generalized 
Mean  Error 

L2  Error 

Norm  (FE) 

CNM-FD 

CNM-FE 

-6.8464  x 10"4 
9.3977  x 10“4 

9.4432  x IO-4 
1.2464  x 10-3 

5.1987  x 10-3 
6.4337  x 10-3 

1.9695  x 10-3 

OIM-FD 

OIM-FB 

1.3598  x 10"4 
1.3598  x 10-4 

1.6546  x 10~4 
1.6546  x 10"4 

1.0025  x 10-3 
1.0025  x 10-3 

1.4031  x IQ"3 

FIM-FD 

FIM-FE 

-5.9604  x 10-3 
-4.1337  x IQ-3 

7.8714  x 10~3 
5.5091  X 10“3 

4.1075  x 10-2 
2.8901  x IQ-2 

2.7448  x IQ-3 

♦See  Table  III 
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Table  VI 


Error  Comparisons  for  the  Various  Methods  for  9 = .08  in  the  Primary  Problem 
and  where  the  Exact  Analytical  Solution  has  been 
substituted  at  the  First  Time  Step. 


Method* 

Pointwise 
Error,  x » .1 

Maximum  Error 
at  Any  Node 

Generalized 
Mean  Error 

L2  Error 

Norm  (FE) 

CNM-FD 

CNM-FE 

-3.0351  x 10‘4 
3.2920  x 10-4 

8.5043  x 10'4 
8.8078  x 10*4 

5.6708  x 10~3 
5.9831  x 10"3 

1.5442  x 10'3 

OIM-FD 

OIM-FE 

1.4103  x lO-5 
1.4103  x 10“5 

2.1640  x lO"5 
2.1640  x lO’5 

1.5280  x 10-4 
1.5280  x 10"4 

9.3405  x 10“4 

FIM-FD 

FIM-FE 

-2.2702  x 10'3 
-1.6021  x IQ’3 

5.8835  x 10'3 
4.2338  x 10' 3 

4.0493  x 10-2 
2.8906  x lO'2 

2.3705  x 10"3 

*See  Table  III 


The  most  striking  feature  of  Tables  III  through  VI  is  that  the 
optimum  implicit  methods  for  both  finite-differences  and  finite- 
elements  give  identical  accuracy.  In  all  four  tables,  the  finite- 
element  method  shows  the  optimum  implicit  method  to  be  more  accurate 
than  the  Crank-Nicolson  or  the  fully  implicit  methods;  while,  in 
the  finite-difference  method.  Tables  III  and  IV  (discontinuous  initial 
conditions)  show  greater  accuracy  for  the  Crank-Nicolson  method,  but 
Tables  V and  VI,  in  which  the  initial  and  boundary  conditions  have 
been  matched,  show  greater  accuracy  for  the  optimum  implicit  method. 
The  Crank-Nicolson  method  appears  to  be  invariably  more  accurate  for 
finite-differences  than  for  linear  finite-elements.  For  the  fully 
implicit  method,  the  finite-element  version  was  the  most  accurate 
in  each  case. 
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The  large  number  of  plots  generated  during  this  investigation 
will  be  analyzed  in  the  order  in  which  they  appear  in  Appendix  H. 
This  appendix  is  divided  into  four  sections.  The  first  deals  with 
the  primary  problem  as  solved  by  finite-differences.  The  second 
deals  with  the  primary  problem  as  solved  by  finite-elements.  The 
third  presents  the  results  of  the  secondary  problem  as  solved  by 
finite-differences.  And  the  fourth  deals  with  the  secondary  problem 
as  solved  by  finite-elements.  The  introductory  remarks  in  the 
appendix  expl&in  the  rationale  for  the  use  of  the  various  error 
norms  and  graphical  formats.  A special  note  is  in  order,  however, 
on  the  interpretation  of  the  discretization  error  ratio.  If  the 
error  for  some  norm  is  given  by 


5(Ax): 


(136) 


then  the  result  of  decreasing  the  interval  size  by  a factor  of  1/2 
will  be 


CjfAx)2 


C 137) 


Similarly,  the  effect  of  halving  the  interval  size  when  the  error  is 
given  by 


5 (Ax)' 


(138) 
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is 


ex  .tjCAx)4 

% " ^(f)4 


16 


(139) 


Thus,  a discretization  error  ratio  of  4 indicates  0(Ax)2  accuracy; 
while  one  of  16  indicates  0(ax)4  accuracy. 

Section  I Results.  This  section  gives  the  results  of  the  primary 
problem  as  solved  by  the  finite-difference  method.  The  first  set  of 
graphs  in  this  sections.  Figures  H-l  through  H-15,  show  the  order  of 
convergence  for  the  case  p = 0.5  when  the  exact  analytical  solution 

has  been  substituted  at  the  first  time  step.  The  first  six  of  these 
show  the  discretization  error  ratio  (DER)  plotted  against  a number  of 
values  of  the  parameter  a . The  peak  occurs  at  or  near  the  value 
of  o predicted  by  Crandall  (a  = .33333)  . Because  of  transient 
effects  on  the  peak  shape  and  position,  the  DER  was  also  plotted  against 
time  in  Figures  H-7  through  H-12  for  three  methods:  optimum  implicit, 
Crank-Nicolson,  and  fully  implicit.  In  each  error  norm,  the  optimum 
implicit  and  the  Crank-Nicolson  methods  approach  0(Ax)2  in  accuracy, 
while  the  optimum  implicit  method  approaches  0(Ax)4  in  accuracy. 

The  behavior  with  respect  to  time  cannot  be  explained.  In  Figures 
H-13  through  H-15,  the  absolute  magnitude  of  the  error  is  seen  to  dip 
to  a minimum  in  the  neighborhood  of  the  predicted  optimum  value  of 
o for  each  error  norm. 
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The  next  three  sets  of  graphs  all  deal  with  results  for  the  case 
p = 1.0  . In  Figures  H-16  through  H-21,  the  numerical  solution  has 

not  been  modified  at  all.  Because  of  the  discontinuity,  the  convergence 
rate  is  not  very  much  enhanced  for  the  optimum  value  of  a (a  = 0.416667)  . 
TTiese  plots  also  show  how  instability  has  begun  to  affect  the  solution 
for  small  values  of  o . In  Figures  H-22  through  H-36,  the  numerical 
solution  has  been  replaced  by  the  exact  analytical  solution,  and 
again,  the  results  approach  the  predicted  behavior,  especially  for  the 
finer  discretization.  In  Figures  H-34  through  H-36,  the  absolute 
magnitude  of  each  error  norm  can  be  seen  to  fall  to  a minimum  in 
the  neighborhood  of  the  predicted  optimum  value.  Figures  H-37  through 
H-42  show  that  the  order  of  convergence  can  be  improved  over  the 
unmodified  solution  by  subdividing  the  space-time  grid  for  the  first 
time  step  as  discussed  previously.  The  improvement  is  not  dramatic, 
however.  More  improvement  could  be  made,  perhaps,  by  using  one  of 
the  other  methods  considered  in  previous  chapters. 

In  the  last  set  of  graphs  in  this  section.  Figures  H-43  through 
H-57,  the  exact  solution  has  been  substituted  at  the  first  time  step 
and  the  peaks  in  the  DER  versus  a curves  are  well  defined  and  fall 
in  the  neighborhood  of  a = 0.45833  , the  predicted  optimum  for 

p = 2.0  . Again,  the  DER  versus  time  curves  show  that  the  DER 

approaches  the  predicted  values  for  each  of  the  methods.  Finally, 

Figures  H-55  through  H-57  indicate  that  the  minimum  absolute  error 
occurs  near  the  predicted  optimum  value  of  a 
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Section  II  Results.  This  section  deals  with  the  primary  problem 


as  solved  by  finite-elements.  The  first  set  of  graphs  in  this 
section  was  generated  for  the  case  p = 0.5  . In  Figures  H-58 

through  H-77,  the  exact  analytical  solution  has  been  substituted  at 
the  first  time  step  to  eliminate  the  problem  with  the  discontinuity. 

Just  as  in  the  finite-difference  method,  each  error  norm  is  seen  to 
approximate  the  predicted  behavior.  This  is  perhaps  best  demonstrated 
in  Figures  H-66  through  H-73.  The  norm  in  Figures  H-64  and  H-65 

is  seen  to  have  a very  broad  peak,  but  in  terms  of  order  of  convergence, 
there  seems  to  be  hardly  any  advantage  for  one  method  over  another. 

Also,  as  in  the  finite-difference  case,  Figures  H-74  through  H-76 
show  that  error  in  the  pointwise  sense  has  a minimum  in  the 
neighborhood  of  the  predicted  optimum  value,  a = 0.66667  . Figure  H-77 

shows  that  the  minimum  in  the  norm  is  a function  of  time  and  occurs 
at  a point  close  to  the  fully  implicit  method. 

The  next  four  sets  of  graphs  all  give  results  for  p = 1.0 
where  the  predicted  optimum  value  of  a is  0.58333.  In  Figures  H-78 
through  H-85,  no  modifications  have  been  made  to  correct  for  the 
effects  of  the  initial  discontinuity.  The  order  of  convergence  seems 
at  times  to  have  a peak  and  at  times  not  to  have  a peak.  Any  peaks 
which  do  appear  are  only  transient  and  occur  for  values  of  o much 
higher  than  the  predicted  optimum  value.  For  Figures  H-86  through  H-105, 
the  exact  analytical  solution  has  been  substituted  at  the  first  time 
step.  As  before  when  this  approach  was  taken,  the  predicted  results 
are  seen.  Also,  in  Figure  H-105,  the  norm  is  seen  to  have  a fairly 
well  defined  minimum,  but  at  a value  of  a which  is  displaced  to  the 
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right  of  the  optimum  value.  In  the  next  set  of  plots.  Figures  H-106 
through  H-113,  an  attempt  was  made  to  reduce  the  effect  of  the  initial 
discontinuity  by  a subdivision  of  the  space-time  grid  for  the  first 
time  step.  Any  improvement  in  the  order  of  convergence  is  only  transient, 
or  at  best,  not  very  significant.  Another  method,  as  suggested  by 
Figures  13  and  14  in  Chapter  III,  is  shown  in  Figures  H-114  through 
H-121.  While  this  method  looked  promising  when  it  was  constructed, 
the  results  were  very  discouraging.  These  figures  show  virtually  no 
advantage  for  any  of  the  stable  methods  in  any  error  norm. 

The  last  set  of  graphs  in  this  section  gives  results  for  the  case 
where  p = 2.0.  In  these  plots.  Figures  H-122  through  H-141,  the  exact 
solution  has  once  again  been  substituted  for  the  numerical  solution  at 
the  first  time  step.  As  before  when  this  was  done,  the  DER  has  its 
peak  at  the  predicted  value  of  a , 0.541667  , and  the  DER  versus 

time  curves  indicate  that,  in  the  norms  which  measure  pointwise  or 
related  error,  the  predicted  convergence  rate  is  approached  for  each 
of  the  three  indicated  methods.  Further,  and  perhaps  even  most 
importantly,  the  L ^ error  norm  in  Figures  H-128  and  H-129  shows  a 
peak  structure  which  is  very  similar  to  that  observed  for  the 
pointwise  error  norms.  This  is  a very  encouraging  development  which 
will  be  more  thoroughly  discussed  in  the  next  chapter.  Finally, 

Figures  H-138  through  H-141  show  that  the  minimum  point  in  each  of  the 
error  norms  is  very  well  defined. 

Section  III  Results.  This  section  deals  with  the  secondary  problem 
as  solved  by  finite-differences , The  analysis  in  this  section  and 
the  next  was  limited  by  time;  therefore,  only  one  value  of  the  Fourier 
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modulus,  p * 1.0  , and  only  two  discretization  schemes,  Ax  = 0.1 

and  Ax  * 0.05  , were  considered. 

In  the  first  set  of  plots,  no  modifications  have  been  made  to  the 
numerical  scheme.  In  this  set.  Figures  H-142  through  H-144,  the 
instabilities  in  the  lower  values  of  a are  apparent,  perhaps  excited 
by  the  discontinuity  between  the  initial  and  boundary  conditions 
since  they  do  not  appear  in  the  case  where  the  exact  solution  is  sub- 
stituted at  the  first  time  step.  In  Figure  H-142,  there  is  a peak 
structure,  but  it  is  displaced  to  the  right  of  the  predicted  optimum 
value  of  a , 0.41667  . 

In  the  second  set  of  plots,  the  exact  analytical  solution  has 
been  substituted  at  the  first  time  step.  This  procedure  shows  the 
sharper  results  which  appear  when  the  discontinuity  is  removed  from 
the  problem.  In  Figures  H-145  through  H-153,  the  same  trends  which 
were  observed  in  the  primary  problem  are  seen  here.  The  minimum 
points  in  Figures  H-152  and  H-153  are  also  very  well  defined. 

In  the  last  set  of  plots.  Figures  H-154  through  H-156,  the 
space-time  grid  has  been  subdivided  for  the  first  time  step,  in  order 
to  try  to  improve  the  convergence  rate  for  the  optimum  method.  Some 
improvement  is  seen  in  the  first  two,  but  in  the  generalized  mean,  which 
computes  the  error  from  nodal  points  spread  over  the  entire  spatial 
domain,  the  improvement  is  seen  to  be  only  slight. 

Section  IV  Results.  This  section  gives  the  results  of  the  secondary 
problem  as  solved  by  the  finite-element  method.  As  in  the  last  section, 
the  analysis  is  restricted  to  the  cases  where  p = 1.0  and  Ax  = 0.1 
and  0.05  . 
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The  first  set  of  plots.  Figures  H-157  through  H-160,  shows  the 


usual  poor  convergence  behavior  of  the  solution  when  there  is  no 
attempt  to  correct  the  problem  with  the  discontinuity  between  the 
initial  and  boundary  conditions.  There  is  a very  poorly  defined 
minimum  in  the  L-  norm  at  around  a * 0.85  . 

In  the  second  set  of  plots.  Figures  H-161  through  H-172,  the  exact 
analytical  solution  has  been  substituted  for  the  first  time  step  and 
predicted  results  once  more  begin  to  appear.  It  is  very  interesting 
to  note  that  the  norm  has  a fairly  well  defined  minimum  for 

o * 0.675 


In  the  last  set  of  plots,  Figures  H-173  through  H-176,  the  space- 
time  grid  was  subdivided  for  the  first  time  step  in  an  attempt  to 
alleviate  the  effect  of  the  initial  discontinuity.  There  is  a nice 
peak  in  the  DER  versus  a curve  for  the  pointwise  error  at  x = 0.1  , 

but  unfortunately  it  occurs  for  a value  of  a which  is  displaced  to 
the  right  of  the  optimum  value.  Otherwise,  there  appears  to  be  no 
great  advantage  in  the  optimum  method  in  this  case,  although  there 
is  some. 

Summary.  A number  of  trends  has  appeared  throughout  the 
foregoing  discussion.  First,  whether  finite-elements  or  finite-differences 
are  used  to  solve  the  problem,  the  discontinuity  between  the  initial 
and  boundary  conditions  causes  problems  with  respect  to  convergence 
for  the  optimum  method.  Second,  the  attempt  to  correct  this  problem 
by  subdividing  the  space-time  grid  for  the  first  time  step  is  helpful, 
but  far  from  satisfactory.  Third,  the  substitution  of  the  exact 
analytical  solution  at  the  first  step  eliminates  the  discontinuity 
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in  the  problem  by  transforming  it  into  a new  problem,  and  for 
this  new  problem  the  difficulties  with  the  convergence  rate  for  the 
optimum  method  disappear.  Fourth,  the  L ^ norm  shows  a minimum  for 
the  finite-element  method  especially  for  the  larger  values  of  the 
Fourier  modulus. 
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Conclusions  and  Recommendations 


Conclusions 


Stability  Analysis.  From  the  results  of  the  stability  analysis, 
it  is  clear  that  the  finite-difference  method  is  more  stable  in 
general  than  the  finite-element  method  with  linear  elements.  It  is 
certainly  curious  that  the  stability  curve  for  the  optimum  finite- 
difference  method  coincides  with  the  one  for  the  optimum  method  in 
linear  finite-elements.  More  will  be  said  about  this  later. 

Ettot  Analysis.  There  are  a number  of  conclusions  and  inferences 


in  this  section.  If  the  problem  considered  has  a discontinuity  of 
the  type  considered  in  the  primary  problem,  then  of  all  the  methods 
of  solution  considered,  the  finite-difference  version  of  the  Crank- 
Nicolson  method  gives  the  most  accurate  results.  This  conclusion 
may  not  hold  true  if  some  acceptable  alternate  handling  of  the 
discontinuity  is  used,  such  as  the  one  suggested  by  Smith  (Ref  8:49). 
For  problems  which  have  no  singularities,  the  optimum  implicit 
method,  by  either  finite-differences  or  finite-elements,  gives  the 
most  accurate  results.  The  fact  that  the  error  is  identical  for 
the  finite-element  version  of  the  optimum  implicit  method  and  the 
finite-difference  version,  or  Crandall's  method,  coupled  with  the  fact 
that  the  stability  curves  for  these  two  methods  coincide,  leads 
to  the  inescapable  conclusion  that  they  are  in  fact  the  same  method, 
only  derived  by  two  widely  different  approaches.  Said  another  way, 
linear  finite-elements  and  finite-differences  coincide  for  the  optimum 
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implicit  method  and,  moreover,  share  an  0(Ax)6  scheme  as  derived 
in  Appendix  D for  finite-elements  and  by  Crandall  for  finite- 
differences  (Ref  2:320), 

In  considering  only  the  Crank-Nicolson  method,  the  finite- 
difference  version  was  more  accurate  in  every  case  than  the  linear 
finite-element  version.  But  it  should  be  borne  in  mind  that  the 
accuracy  of  the  finite-element  version  can  be  enhanced  not  only 
by  further  refinement  of  the  space  mesh,  but  also  by  increasing  the 
order  of  the  polynomials  used  for  the  approximate  temperature 
distributions  within  the  elements.  The  previous  conclusion,  then, 
is  largely  qualified  by  the  restriction  to  linear  elements.  However, 
since  the  time  derivatives  are  approximated  by  a finite-difference 
expression  for  which  no  parallel  increase  in  accuracy  can  be  attained 
by  going  to  higher-order  polynomial  approximate  temperature  distributions, 
then  there  may  be  a serious  restriction  on  the  total  increase  in 
accuracy  which  can  be  attained  in  this  way.  This  problem  may  be 
solved,  however,  by  using  a variational  principle  such  as  the  one 
proposed  by  Curtin  (Ref  11:255). 

It  is  suggested  that  the  problems  with  convergence  for  the 
optimum  schemes  are  due,  at  least  in  part,  to  the  weight  given  to  the 
higher  eigenfunctions  in  the  expansion  of  the  initial  condition  in 
the  case  of  a discontinuity.  These  higher  eigenfunctions  are 
automatically  filtered  out  and  lost,  since  the  Crank-Nicolson  method 
is  nothing  more  than  a Padd  rational  approximation  for  the  temporal 
behavior  of  the  solution  (Ref  17:266): 
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exp  (-XA0) 


(140) 


1 - X(A8/2) 
1 ♦ X (A0/2) 


thus  giving  rise  to  difficulty  with  problems  for  which  these  higher 
eigenfunctions  are  important.  If,  however,  a more  accurate  Padd 
rational  approximation  is  used,  for  example 


exp  (-XA0) 


1 - X (A8/2)  ♦ X2(A82/12) 

1 ♦ X (A0/2)  ♦ X2 (A02/12) 


(141) 


then  perhaps  fewer  of  the  higher  eigenfunctions  will  be  lost  and  can 
contribute  to  the  numerical  solution. 

It  is  very  interesting  to  surmise  on  the  possibility  of 
attaining  a more  accurate  solution  for  finite-elements  in  the 
norm.  The  results  indicate  a trend  that,  as  the  Fourier  modulus  p 
gets  larger,  the  DER  versus  a curves  for  the  L2  norm  look  more 
and  more  like  those  for  the  pointwise  norms.  This  would  mean  that  there 
is  a minimum  error  condition  for  large  p in  the  norm  which  is 

perhaps  predictable.  This  would  be  of  great  benefit,  since  an  engineer 
using  such  a scheme  would  have  an  approximate  solution  at  every  point 
on  the  domain,  not  just  at  the  nodes,  and  know  with  some  degree  of 
confidence  that  the  solution  is  accurate  to  a high  degree  everywhere 
on  the  domain. 

From  the  numerous  plots  in  Appendix  H,  it  is  clearly  shown  that, 
for  a problem  without  singularities,  the  predicted  rates  of  convergence 
are  approached  and  the  minimum  error  condition  occurs  for  the  predicted 
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I 

I 

I 


values  of  o . This  is  true  for  both  the  finite-difference  and  finite 
element  formulations. 

Recommendations 

As  with  many  investigations,  this  project  uncovered  more 
questions  than  it  answered.  For  one,  it  would  be  interesting  to  see 
the  result  of  using  a higher-order  polynomial,  or  perhaps  splines, 
for  the  elemental  temperature  distribution  coupled  with  the  use  of  .1 
more  accurate  Fadd  rational  approximation  for  the  temporal  behavior  of 
the  solution.  Varga  (Ref  18:70)  suggests  the  use  of  Tchebycheff 
semidiscrete  approximations  to  spread  the  error  in  the  approximation 
of  exp  (-XA0)  over  the  entire  time  domain.  It  would  be  interesting 
to  investigate  these  modifications.  Additionally,  further  work 
could  be  done  in  two  and  three  dimensions  to  try  to  extend  the  method. 
It  would  also  be  of  interest  to  investigate  the  error  at  the  right 
boundary  for  the  secondary  problem,  since  a lower  order  of  accuracy 
is  predicted  at  that  point.  Last,  and  probably  most  important,  the 
behavior  of  the  norm  for  large  values  of  the  Fourier  modulus, 

p , should  be  carefully  investigated,  as  this  could  lead  to  an 
optimum  method  in  terms  of  this  norm,  one  which  is  fundamental  to 
the  finite-element  concept. 
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Appendix  A 


The  Analytical  Solution  for  the  Primary  Problem 


In  normalized  form,  the  primary  problem  in  this  investigation  was 


3u  3fu 
30  " ...2 


(A—  1) 


u(x,0)  ■ 1,  0*0 


(A-2) 


u(O,0)  * u(l,0)  =0  , 0 > 0 


(A-3) 


where 


* normalized  tempera tue 

* normalized  time  variable 

* normalized  spatial  variable  (the  bar  has 
been  dropped  from  x for  simplicity) 


The  solution  will  be  developed,  using  the  method  of  separation 
of  variables.  Assume  u(x,0)  * X(x)0(0)  . Then 


•i  * 


(*-<) 


8 X- 


(A-5) 


3U  y dO  A y 

5F  - * 3?  " ‘ 9 


Substitution  of  these  results  into  Equations  (A-l)  through  (A-3)  gives 


(A-6) 


(A-7) 


since  X is  a function  of  x only  and  G is  a function  of  0 only, 
and  where  y is  a separation  constant.  Separation  of  the  boundary 
conditions  gives 


X(0)  9 (61 


6 > 0 


(A-8) 


x(i)  e (0) 


0 > 0 


(A-9) 


Since  Equations  (A-$l  and  ( A-9'l  must  hold  for  any  0 > 0 , then 


(A-10) 


With  the  boundary  conditions  of  Equation  (A- 101  the  Sturm-liouville 
problem  contained  within  (A-l)  through  (A-3)  is  given  by 


X"  - yX 


(A- 11) 
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X(0) 


x(i) 


(A-12) 


The  solution  of  the  Sturm-Liouville  problem  must  be  broken  into 
three  cases.  In  the  first  two  cases,  where  y > 0 and  y = 0 
respectively,  there  is  no  solution.  In  the  third  case,  y * 0 , let 
Y * -a2  where  a f 0 , then  solving  the  characteristic  equation  for 

the  differential  equation  (A- 11)  gives 

X(x)  = A cos  (ax)  ♦ B sin(ax) 

as  the  general  solution.  The  boundary  conditions  are  then  substituted 
to  obtain 

X(0)  - A • 1 + B • 0 = 0 


or 


and 


X(l) 


(0)  . 1 + B sin(c) 


(A-13) 


(A-14) 


(A-15) 


(A- 16) 


or 


sin(a)  * 0 


(A-17) 
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since  B cannot  be  zero  if  a non-trivial  solution  is  desired.  There 


< 


k 


( 


is  an  infinite  number  of  values  of  o which  will  satisfy  Equation  (A- 17) 
which,  after  it  is  recalled  that  a / 0 and  that  only  positive  values  are 
needed  since  the  negative  values  only  repeat  the  solutions  given  by  the 
positive  values,  are  given  by 

- mr  , n = 1,2,3,  •••  (A-18) 

Now  the  eigenvalues  of  Equation  (A- 11)  are  given  by 

Y„  = -(on)2  = -(mr)2  , n = 1,2,3,...  (A-19) 

The  solutions  of  Equation  (A- 13)  are  then 

Xn(x)  = Bn  sin(mrx)  , n = 1,2,3,...  (A-20) 

In  the  second  part  of  the  problem 

0'  - Y0  =0  (A-21) 

After  integrating.  Equation  (A-21)  becomes 

0(0)  = Cn  exp(-(mr)20)  (A-22) 

By  invoking  the  principle  of  superposition,  the  general  solution  of 
Equation  (A-l)  is  given  by 


& 
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(A-23) 


u(x,6)  - l (B  * C ) sin(nirx)  exp(-(nn)20) 

n-1  n 


The  initial  conditions  are  then  expanded  in  an  infinite  Fourier  series 
of  the  eigenfunctions.  Equation  (A-20) , as  follows: 


u(x,0)  > [ A sin(nnx)  *1  * 1 

n-1 


where  A * B • C , and  zero  has  been  substituted  for  theta, 
n n n 

This  should  be  recognized  as  a Fourier  sine  series  where 


(A—  24) 


\ * 


(1  * Xn(x)) 

TiyiiF 


(A-25) 


in  inner  product  notation  introduced  earlier  and 


|Xn(x)|'2  " / XS(x)  dx 


or  since 


(1  • Xn(x)) 


(A-26) 


(A-27) 
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and 


( 


l|XnCx}||2  - i 


(A-28) 


then 


A„  ■ ST  - (-»") 


(A-29) 


If  n » 2m  - 1 , then  Equation  (A-23)  can  be  written 


u(x,0)  * l 4 sin((2m-l)irx)  exp(-((2m-l)ir)20)  (A-30) 

m-1  (2m-l)ir 


which  is  the  solution. 
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Appendix  B 

The  Analytical  Solution  for  the  Secondary  Problem 


In  normalized  form,  the  secondary  problem  is  given  by 


3u 

30 


lln. 

3x2 


(B-l) 


u(o,e) 


o,  e > o 


(B— 2) 


u(x,6)  * 1 , 0=0 


(B-3) 


As  in  Appendix  A,  the  method  of  separation  of  variables  will  be  used. 
A product  type  solution,  u(x,0)  = X(x)  0(0)  will  be  assumed.  After 
substitution  into  Equations  (B-l)  through  (B-3),  the  problem  becomes 


xr^xi 

X(x) 


0*(9) 

0(0) 


* Y 


(B-4) 


where  the  prime  indicates  differentiation  with  respect  to  the 
indicated  independent  variable,  and  y is  a separation  constant. 
Similarly,  the  boundary  conditions  become 

X(0)  = X'(l)  = 0 (B-5) 


104 


The  Sturm-Liouville  problem  associated  with  Equation  (B-l)  is  then 


X~  - yX  * 0 


X(0)  - X'(l)  ■ 0 


(B-6) 


CB-7) 


For  the  cases  where  y > 0 and  y = 0 , there  is  no  solution  to 

Equations  (B-6)  and  (B-7) . In  the  last  case,  where  y < 0 , if 


y * -a2  , a i 0 


(B-8) 


then  the  solution  of  Equation  (B-6)  is  given  by  solving  the  resulting 
characteristic  equation  to  give 


X(x)  = A cos (ox)  + B sin (ax) 


(B-9) 


After  substitution  of  the  boundary  conditions,  the  result  is 


X(0)  = A • 1 + B * 0 


(B-10) 


A ■ 0 


(B-ll) 
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and 


X'(l)  ■ o B cos(a)  * 0 

or,  since  a f 0 and  B cannot  equal  zero  if  the  solution  is  to 
be  non-trivial 

cos (a)  = 0 

The  values  of  a which  satisfy  Equation  (B-13)  are 


nj  , n * 1,3,5, 


or  if  n » (2m- 1) 


am  * (2m-l)-  , m = 1,2,3, 


(B-12) 


(B-13) 


(B- 14) 


(B— 15) 


The  eigenvalues  of  the  problem  are  given  by 


Ym  * -(ara)2  * -((2m-l)J)2  , m a 1,2,3,.,,  (B-16) 


The  solutions  of  Equation  (B-6)  with  boundary  conditions  (B-7)  are  then 


- Bm  sin((2m-l)J  x) 


(Bt17) 
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In  the  second  part  of  the  problem 


e'(e)  - y 0 


integration  with  respect  to  0 yields 


Qmce)  - Cm  exp(-C(2m-l)j)2  0) 


The  use  of  the  principle  of  superposition  then  leads  to 


u(x,e)  * l Affl  sin((2m-l)j  x)  exp(-((2m-l)j)28) 

n*l 


where  A ■ B «C 


m mm 

Application  of  the  initial  conditions  gives 

00 

u(x,0)  = l Am  sin((2m-l)J  x)  • 1 = 1 

m*l 

The  coefficients,  A , are  recognized  as  the  coefficients  of  the 

in 

Fourier  sine  series  given  by 


A» 


l|Xm(x)||2 


(B- 18) 


(B- 19) 


(B-20) 


(B-21) 


(B-22) 
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where 


1 


* x_ 


(x)) 


/ sin  ((2m-l) j x)  dx 
o * 


(2n-l)it 


(B-23) 


and 

1 . 

HX^CxJll2  > / sin2((2m-l)j  x)  dx  - j (B-24) 

Therefore,  the  solution  of  Equation  (B-l)  with  boundary  conditions  (B-2) 
and  initial  condition  (B-3)  is 

m 

u(x,B)  - y 4 sin((2m-l)T  x)  exp(-((2m-l)T)20)  (B-25) 

»-l  "'■(HPTTir  2 2 
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Appendix  C 


a i i mp 


Derivation  of  the  Truncation  Error 
for  the  Finite-Difference  Formulation 


Crandall  (Ref  2:319)  outlines  a method  for  investigating  the 
truncation  error  of  Equation  (28).  The  details  are  given  by  Campbell 
(Ref  15:89),  with  a few  minor  errors,  for  the  case  which  includes  a 
diffusivity  constant.  Here,  the  correct  derivation  will  be  given  for 
the  normalised  equation  which,  of  course,  carries  the  diffusivity 
constant  implicitly. 

In  Equation  (28),  the  exact  solution  u(x,6)  has  been  expanded 
ov.r  six  points  . »l>kM  . V,  >t.,  . u,  , k . 

u^  ^ , and  u^  ^ , The  difference  between  Equation  (28)  and  the 
exact  differential  equation  is  called  truncation  error  and  is  given  by 


(Vl,k  + l ~ ^i.k+l 


. (uitk,+i  -._v0 

A0  Ax2 


u.  , . - 2u,  ♦ uSj,  . 

( 1 -a)  ( * i»k  i*l.k^ 


♦ u 


i •*■  1 , k ♦ 1 ) 


Ax* 


3u  02u 

(36  ' 0X2 


(C-l) 


Since  3u  a2u  n , then 
30  'ax2 


t * &a  lAi  ui«l,k*l  + A2  ui,k*l 


♦ A.  u 


1 “i*l,kel 


B1  Ui-l,k  ” ®2  Ui,k  " Bl  “i+l.k1 


(C-2) 


1 


I 
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where 


1 ♦ 2 pa 


P(l-a) 


1 - 2p(l-a) 


It  is  a standard  shorthand  convention  to  let  Ax  = h and  A6  = k 
If  this  notation  is  used,  the  six  points  of  u , mentioned  before, 
can  be  expanded  in  a Taylor  series  centered  about  u.  . to  give 

1 »K 


'i.k*! 


"i.k  * kue  * T U20  * lru39  *■ 


(C-3) 


"Ul.k  ■ "i.k  1 hux  * Tu2x  1 JT  uJx  *••• 


(C-4) 


and,  if  k • ph2  , then 


u.xl  = u.  , ± hu  + h2(p  i-)  u 

i±l,k+l  i,k  x 2 2x 


± l>3  (P  ♦ 5)  u3x  * b1*  4 * f * & U4x 


‘ kS  (V  * 6 * ®>  u5x 


(C-5) 


where  the  subscripts  indicate  differentiation  with  respect  to  the 
indicated  independent  variable  for  the  indicated  number  of  times  at  the 
point  (iAx,kA0)  . If  Equations  (C-3)  through  (C-5)  are  substituted 
into  Equation  (C-2)  and  if  the  following  relations  are  used 


u6  = u2x 


(C-6) 


U20  = U02x  = U4x 


(C-7) 


U30  = U202x  = U04x  = U6x 


(C-8) 


along  with  the  fact  that  k = ph2 « then  each  power  of  h may  be 
collected  together  in  the  truncation  error.  The  coefficients  of  all 
of  the  odd  powers  of  h cancel  to  zero.  The  coefficients  of  the  even 
powers  of  h also  cancel  to  zero  up  to  and  including  order  two.  The 
resulting  expression  is 
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or 


*t  ' »2  (f  - ap  - i)  u4x 


* h'4  -^-5  - W>“6X  * 


The  truncation  error  is  then  of  order  h2  unless 


(f  - «P  - tj) 


or 


T (1 ' 


(C-10) 


(C-11) 


CC— 12) 


This  last  expression  is  the  one  derived  by  Crandall.  Crandall  also 
points  out  that,  if  Equation  (C-12)  is  substituted  into  the  coefficient 
of  the  fourth  power  of  h term,  and  the  resulting  expression  set  equal  to 
zero,  then  the  result  would  be  an  order  h6  expression.  The  values 
of  a and  p which  yield  this  0(h)6  scheme  are  (Ref  2:320) 


£ 

10 
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CC*13) 


and 


P 


3-/5 

6 


(C-14) 


which  is  merely  the  point  of  intersection  of  the  coefficients  of  the 
h2  and  h4  terms  in  the  truncation  error. 

Dirichlet  boundary  conditions  do  not  alter  this  conclusion  since, 
for  the  boundary  points,  the  Taylor  series  expansions  are  given  by 


u. 


i-l,k 


u. 


i,k 


hu  + 
x 


h2 

2 


U2x  +' 


= 0 


(C-15) 


and 


u. 


i-l,k+l 


u.  i. 
i»k 


- hu 


h2  (p 


* 2>  U2x 


- 0 (( 


So,  although  the  exact  value  of  the  function  at  the  points  ((i-l)Ax,kA9) 
and  ((i-1) Ax, (k+l)A9)  are  known,  this  does  not  affect  the  Taylor 
expansions  or  the  result  of  their  use  in  Equation  (C-2) . 


-16) 


Appendix  D 


Derivation  of  the  Truncation  Error 
for  the  Finite-Element  Formulation 


The  same  method  which  was  used  in  Appendix  C will  be  used  for  the 
truncation  error  in  the  finite-element  formulation. 

If  the  system  of  equations  (106),  where  u^  = 0 , is  treated 

as  if  it  were  merely  a set  of  difference  equations,  then  the  general 
expression  is 


! 


A,  u.  , , . + A„  u.  . , + A,  u.  , , , = B,  u.  , . + B„  u.  , + B,  u.  , . 

1 i-I, k+1  2 i,k+l  1 l+l, k+1  1 l-l, k 2 i,k  1 i+l,l 


where 

Aj  = 1 = 6ap 

A2  = 4 + 12ap 

Bj  = 1 + 6(l-a)p 

B2  = 4 - 12 (l-a)p 

The  truncation  error  is  given  by  the  difference  between  the  exact  partial 
differential  equation  and  the  difference  equation  (D-l) : 
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t A0 


^A1  Ui-l,k+l  + A2  Ui,k+1  + A1  ui+l,k+l 


- B1  ui-l,k  ' B2  “i,k  - B1  “i.l.k  * <1?  - (D-2> 

oX 


But  since 


. !iu  * o 

30  3x2 


(D-3) 


then 


i = ——  [A,  u.  , , , + A_  u.  , , ♦ A,  u.  , 

t A0  1 1 l-l, k+1  2 i,k+l  1 i+l. 


k+1 


‘ BlVl,k  - B2Ui,k  - BlUi+l,kl 


CD-4) 


Now  if 


ph2 


(D-5) 


U0  = U2x 


CD— 6) 


U20  = U02x  = U4x 


(D-7) 
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(D-8) 


30 


20  2x 


04x 


6x 


where  k = A0  , h ■ Ax  and  the  subscripts  indicate  differentiation 
with  respect  to  the  subscript  and  for  the  indicated  number  of  times  at 
the  point  (iAx,kA0)  , then  the  Taylor  expansions  in  Equations  (C-3) 
through  (C-5)  can  be  substituted  into  Equation  (D-4)  tc  give 


>t  ■ »a  - »P  ♦ T2>  “4x 


♦ h4  (p2(l-3o)  * p(i  - j)  » -i-)  u6x  . °(h)‘  ....  (D-9) 


To  obtain  an  order  h1*  expression,  let 


\ 


“ ap  + 12  = 0 


(D-10) 


or 


7 (1  * ^ 


CD-11) 


This  is  the  finite-element  equivalent  of  the  method  derived  by  Crandall 

v 

for  finite-differences.  It  should  be  noted  further  that  substitution  of 
Equation  (D- 11)  into  the  coefficient  of  the  h4  term  in  Equation  (P-9) 
and  setting  the  resulting  expression  equal  to  zero,  gives  an  order 
h6  expression: 
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(D-12) 


3 ♦ /S' 


(D- 13) 


The  similarities  between  Equations  (D-ll),  (D-12),  and  (D-13)  and  those 
given  by  Crandall  for  the  order  h4  and  order  h6  schemes  are 
striking.  Also,  as  in  the  finite-difference  case,  Dirichlet  boundary 
conditions  do  not  affect  the  truncation  error. 


Appendix  E 

Derivation  of  the  Variational  Statement 

First,  it  will  be  demonstrated  that  the  approach  suggested  by 
Meyer  (Ref  3:399)  will  lead  to  the  wrong  conclusion  concerning  the 
variational  principle. 

Meyer  suggests  that  the  functional  to  be  minimized  is  of  the  form 

1(G)  - / 1 F (x,0,a,f£,|£)  dx  (E-l) 


where  ii  = u(x,0)  ♦ £v(x,0)  , and  v is  an  arbitrary  function 

which  satisfies  the  essential  boundary  conditions  on  the  problem. 

If  Equation  (E-l)  is  differentiated  with  respect  to  i then,  by  the 
chain  rule  of  the  calculus 


ii  <!£  3F  —*  it-  i a 

H * V3u  H * 3uv  3£  * 3u  H 1 x 

O AO 


Now,  since 


(E-2) 


3Q 

H 


v (x , 0) 


(E-3) 


3v 

ax 


(E-4) 
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and 


3ii0 

TC 


3v 

30 


(E-S) 


then,  after  substitution  of  these  back  into  Equation  (E-2) , the  result 
is 


3u  3x 

X 


3F  3v 

3u.  30 
8 


] dx  (E-6) 


Integration  of  Equation  (E-6) 


by  parts  leads  to 


il 


3F  3v 

30.  30 

8 


dx 

X 


(E-7) 


For  the  expression  in  Equation  (E-7)  to  be  equal  to  zero,  the  expression 
in  the  brackets  must  be  equal  to  zero,  or  since  for  the  minimum 
condition,  £ = 0 , then  u = u and 


v 


d 

d7 


(—  )) 
l3u 

x 


3v  3F 

30  3u„ 

0 


0 


(E-8) 


But  this  expression  must  hold  for  any  arbitrary  function  v(x,e) 
which  satisfies  the  essential  boundary  conditions  imposed  on  the 
problem.  Since  the  function  v and  its  derivative  with  respect  to 
Q cannot  be  removed  from  the  expression,  it  is  worthless  as  a tool 

to  obtain  the  functional  F needed  for  the  variational  statement. 
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The  solution  to  this  problem  is  a simple  one.  At  any  given 
point  in  time,  the  problem  can  be  considered  to  be  elliptic,  thus 
removing  0 and  derivatives  with  respect  to  0 from  the  functional 
F . After  this  is  done.  Equation  (E-l)  becomes 


ICS) 


J V (x,  u,  dx 


(E-9) 


where  u(x)  » u(x)  + £v(x)  for  some  fixed  point  in  time.  If 
Equation  (E-9)  is  differentiated  with  respect  to  £ the  chain  rule  of 
the  calculus  leads  to 


3_I 

H 


1 

/ 


o 


(E- 10) 


And  since  now 


3u 


v(x) 


and 


then,  after  substitution  and  integration  by  parts 


(E- 11) 


(E-12) 
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At  the  minimum  point,  u = u and  £ = 0 and  also 


II  - 0 

35 


CE— 14) 


For  Equation  (E-14)  to  hold,  the  expression  in  the  brackets  must  hold 
for  any  arbitrary  v(x)  which  satisfies  the  essential  boundary 
conditions.  Thus 


3F  3 ,3F  . . 

an  " -57  (^ir)  * 0 

x 


(E-15) 


This  equation  is  known  as  the  Euler-Lagrange  equation.  A comparison 
of  Equation  (E-15)  and  the  differential  equation 


(E— 16) 


31  t r3F  , . , . 3 ,3F 

H ■ > [5a  vW  ' v(x)  57  ix 


(B-13) 


which  can  be  rearranged  as 


- A 

3x  l3xJ 


leads  to  the  observation  that 


CE— 17) 


(E-18) 
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and 


3F  m 3u 
3u  3x 

X 


Integration  of  Equation  (E- 18)  gives 


F 


+ f(ux} 


(E-19) 


(E-20) 


and  integration  of 


Equation  (E-19)  leads  to 


F 


g(u) 


(E— 21) 


The  functions  f and  g can  be  found  by  comparing  Equations  (E-20) 
and  (E-21)  so  that 

f ■ T t]Rf  <“2’  * (|?)2  1 CE-22) 

The  variational  statement  of  the  differential  equation  is  then 

1 * I / (u2)  + (|l)2  1 dx  (E-23) 

o 

i 

which  is,  after  normalization,  the  same  variational  statement  given  by- 
Meyer  . 
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Appendix  F 

Development  of  the  Method  of  Weighted  Residuals 


The  method  of  weighted  residuals,  or  the  Galerkin  method,  will 
be  developed  by  following  the  treatment  given  a two-dimensional  problem 
by  Yalamanchili  and  Chu  (Ref  16:10-17). 

The  error  incurred  by  using  the  difference  equation 

♦ (1-a)  — (F-l) 

3x2 

i,k+l  i,k 

where  the  subscript  i indicates  the  spatial  node,  and  the  subscript 
k indicates  the  time  step,  to  approximate  the  following  differential 
equation 


i,k+l 


- u 


A0 


Ul 


= a- 


32u 
i 

3x2 


3u 

36 


3x2 


(F-2) 


is  called  a residual,  defined  by 


R(x)  - a— 
3x2 


♦ (1-a) 


3fu 

3x2 


i,k*l 


ui  k+1  ui  k 
- "T e — ^ 


i.k 


The  object  in  the  method  of  weighted  residuals  is  to  minimize  and 
distribute  this  error  over  the  interval  (0,1)  , or 
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/ R(x)  W(x)  dx  =■  0 

o 


(F-4) 


First,  the  neighborhood  of  the  ith  nodal  point  must  be  divided 
into  two  adjacent  elements.  Then,  for  this  case,  a linear  temperature 
distribution  is  assumed  within  each  element.  For  element  one  at 
time  0 * kA0 


u(x) 


u. 


X.  - X 

f 1 


i-l,k  - x 


) + V c 


X - X 


i-1 


i.k  Xi  - Xj.j' 


(F-5) 


and  for  element  two 


u(x) 


X.  , - X 

(-tii — ) 

'V  _ V ' 


i Ir  'y 

l.K  Xl+1 


X - Xi 

ui*l,k  (; 


Xi*l  ' xi 


(F-6) 


The  weights  W(x)  are  chosen  to  be  the  coefficients  of  the  discrete 
temperatures  in  the  element  equations  (F-5)  and  (F-6) . Outside  of 
these  two  one-dimensional  elements,  the  weighting  functions  are  zero. 
If  these  weighting  functions  are  substituted  into  Equation  (F-4)  along 
with  the  residual,  Equation  (F-3) , then 


f 1 r 3*u 
i*l  3 


♦ / l*}a2u 

x,  (V 


♦ (1-cO 


3fu 

3x2 


k*l 


u. 


- (- 


i.k*l 


i k 

A0  =-*-—) ) (; 


X - X 


i-1 


x « x 

i xi-l 


(1-a) — - 
3x2 


k+1 


u.  , . , - u.  . x...  - x 

( 1,k  -rr ^))  C-  T-v  ) (p-7) 


A0 


VX.  , - X.‘ 
1+1  1 
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Next,  the  temperature  distributions.  Equations  (F-5)  and  (F-6),are 
substituted  into  Equation  (F-7)  and  the  following  finite-difference 
expression  for  the  Laplacian 


Vi  • 2ui  * Vi 


(F-8) 


is  also  substituted  into  Equation  (F-7),  then  Equation  (F-7)  becomes 


X • x — X 

^ ~2  ^Ui-l,k+l  ” 2ui,k-*l  * Ui*l,k+1^  ^ Ax"  ^ 

i-1 


X.  - . X - X.  , 

* / ^ (ui-i,k  - 2ui,k ♦ Vi.k)  (-at11) 


xi  1 Ax 


-0  -do tV,.k*l  * “i,k.l 

i-1 


X.  - X X - X.  . 

- ui-i,k  - ui,k  1 dx 


I 1 777  (Ui-l,k+l  " 2ui,k+l  + Ui*l,k+1^  ( 


x.  , - x 
i+l 


Xi  Ax2 


X,  . , - X 


Xj  Ax2 


X“1  mi  <Vi  - ’>  h.k.i  (Hr1)  * V.,W  'Hr1' 


,xi*l  ~ x 
Ui,k  Ax 


x - x. 


) - vi.k  <Tr>  i d*  - 0 


(F— 9) 
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After  integration  and  simplification.  Equation  (F-9)  becomes 


a / 

2Ax  lui-l,k+l 


2ui,k+l  + ui+l,k+l^ 


+ 2A?^  (ui-l,k  “ 2ui,k  * “i+l.k5 


u 

6A0  i-l,k+l 


Ax  Ax  Ax 

3A0  Ui,k+1  + 6A0  Ui-l,k  3A6  Ui,k 


+ 2 Ax  (ui-l,k+l 


2u . . , + u . , . , ) 
i,k+l  l+l.k+l-' 


+ (ui-l,k  " 2ui,k+l  + ui+l,k5 


(F-10) 


Ax  Ax  Ax  Ax  „ 

6A0  i,k+l  ■ 6A0  ui+l,k+l  3A0  i,k  6A0  ui+l,k 


If  p * A0/(Ax)2  , then  Equation  (F-10)  can  be  rewritten  as 


- °P)  “t-l.k.l  * (!  - 2“P>  ui,k  * (i  - ui+l,k+l 

(F-U) 

• (|  * d-a)p)  Vl  k . (f  - 2(l-a)p)  u1>k  • (i  ♦ »-«)P)  Vl>k 

This  is  the  general  finite-element  difference  equation  which,  on  close 
inspection,  is  seen  to  be  exactly  equivalent  to  the  general  expression 
which  was  derived  by  the  Ritz  method  using  the  variational  principle 
given  by  Meyers. 


Appendix  G 

Derivation  of  the  L 2 Error  Norm  for  Finite-Elements 

Strang  and  Fix  (Ref  9:49)  show,  using  the  famous  "Nitsche  trick," 

* 

that  a finite-element  approximation,  with  piecewise  linear  elemental 
temperature  distribution,  satisfies 


| | u - uN| |q  < c/h2  (G-l) 

where  the  norm  is  the  H°  or  norm  defined  in  Chapter  II. 

This  error  measure  is  more  closely  related  to  the  underlying  nature 
of  the  finite-element  method  than  is  a pointwise  error  bound;  it  is 
often  referred  to  as  the  displacement  error. 

The  exact  analytical  solution,  u , is  given  by 


u(x,0)  * l 


n*l  t2"-1)* 


sin((2n-l)irx)  exp(-((2n-l)iv)20) 


(G-2) 


The  finite-element  solution  assumes  a linear  temperature  distribution 
within  the  elements  defined  over  the  interval  (iAx, (i+1) Ax)  : 


N,  . 
u (x) 


0 » kA0 


Vk  * 


-)  (x  - iAx) 


(G-3) 


Substitution  of  Equations  (G-2)  mid  (G-3)  into  the  L2  norm  gives 
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lu-uN|l0  = [ / (u(x)  - UN(X))2  dx  l*5 

O 


N-l  (i+l)Ax  00 

[l  / ([  l a sin(b  x)  exp(-b20)]  - (C  + DJ 

i=0  iAx  n=l  n n n x 


where 


a„  = 4/ (2n-l)  it 


C 

D 

N 


(2n-l)n 

ui,k  - <iAX<ui+i.k  - Ui,k)/Ax) 

(Vi.k  - Ui,k)/Ax 

number  of  elements 


Since  the  integration  is  taken  with  respect  to  x , then 


En  = exp(-b2e) 


as  a further  simplification.  Expanding  Equation  (G-4) 


2 dx  ]**  (G-4) 


let 

(G-S) 
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N N-l  (i*l)Ax  «> 

I lu  - u ll0  * fl  (/  d a sin(b  x)  E )2  dx 

i-0  iAx  n=l  n n n 


(i+ljAx  « 

2/  ( l an  sin(bnx)  EJ  (C  + D^)  dx 


iAx  n=l 


(i+1)  * l, 

♦ / (C  ♦ Dx)2  dx  ) ] 2 

i x 


(G-6) 


This  last  step  can  be  taken  since  it  is  well  known  that  in  the 
following  equation 


‘Is"1  (JiV 


(G-7) 


if  the  series  whose  terms  are  s and  the  series  whose  terms  are  t 

n n 

both  converge  absolutely,  then  the  series  whose  terms  are  wn  also 
converges  absolutely.  The  exact  solution  has  previously  been  shown 
to  be  uniformly  convergent.  Further,  the  integrals  which  involve  infinite 
series  in  Equation  (G-6)  are  integrable  since,  in  the  equation 


S(x)  * l v (x)  , a < x < b 

n=l 


(G-8) 


if  the  functions  vn00  are  continuous  (they  are),  and  if  the  series 
in  Equation  (G-8)  whose  terms  are  vfi(x)  is  uniformly  convergent 
(previously  demonstrated  for  the  first  integral  on  the  right-hand  side 


of  Equation  (G-6),  then  S(x)  is  also  continuous  and  can  be  integrated 
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term  by  term  (Ref  5:28).  The  series 


P(x)  - [ l (an  sinfb  x)  E ) ] (C  + Dx)  = £ Pfx)  (G-9) 

n*l  n n " n=>l 

can  be  shown  to  be  uniformly  convergent  by  the  "Weierstrass  M-test" 

(Ref  5:28).  If  there  is  a convergent  series  of  positive  constants 


T M (G-10) 

L.  " 
n=l 


such  that 


|Pn(x)|  s Mn  , a s x * b 


then  series  (G-10)  is  uniformly  convergent  on  (a,b)  . The  terms 

defined  by  Equation  (G-9)  satisfy  the  relation 


4 

(2n-l)ir 


exp(-((2n-l)7t)20) 


(G-ll) 


since 


uN (x)  < 1 


(G-12) 


and 

|stn((2n-l)irx)  | < 1 


CG-13) 
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Then,  by  the  ratio  test  for  convergence  of  an  infinite  series 


f2n*nV  exP(-(2n+l)2iT29) 
J2i^  exp(-(2n-l)2ir2e) 


0 


(G- 14) 


which  implies  convergence  and,  therefore,  uniform  convergence  for 
Equation  (G-9) . 

Now  the  first  integral  in  Equation  (G-6)  can  be  written 


(i*l)Ax  “ (i*l)Ax  «® 

I,  ■ / (I  an  sin(b  x)  E )2  d x - / l f2(x)  dx  (G-15) 

1 iAx  n-1  n n n iAx  n-1  n 


This  series,  previously  shown  to  be  square-integrable,  can  be  written 


(i*l)Ax 

iAx 


GO 

I f£(x)  dx 
n-1 


» (i*l)Ax 

a I / 

n-1  iAx 


f 2 (x)  dx 
n 


• “ (i*l)Ax 

♦ 2 l [ l J f (x)  f (x)  dx  ] (G-16) 

n-1  m-n+1  iAx 

The  series  is  simply  truncated  when  the  computer  reaches  its  limit  for 
storing  small  numbers.  No  effort  was  made  to  estimate  the  error,  but  it 
could  be  done  by  a method  analogous  to  the  one  used  in  Chapter  II  for 
the  analytical  solution.  After  integration,  the  expression  in 
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Equation  (G-16)  becomes 


ri  5 £ a£E£  [^T  • W-  (sin  2bn(i+1)Ax  - sin  2bn  iA*)  ] 

n»l  n 


* j!  JL  2\w»^(sin^''-b")(inux) 


- Sin((bn-bm)i4x))  - ‘ - (sin((b  *b  ) (i*l)Ax) 


sin((b  +b  )iAx))  ] 
n m 


(G-17) 


where  K is  large. 

The  second  integral  in  Equation  (G-6)  can  be  shown  to  be  termwise 
integrable,  using  the  same  rationale  as  used  for  the  first  integral. 
The  result,  then,  after  integration  is 


K a E 

I2  * 2C  \ -jj — (cos (bn(i+l) Ax)  - cos(b  iAx)) 
n«l  n 


K a E 

-2D  J -Py  — [ (sin(b  (i+l)Ax)  - sin(b  iAx)) 


n-1  n 


- (bn(i+l)Ax  cos(bn(i+l)Ax)  - bniAx  cos(bniAx))  ] 


(G-18) 
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The  third  integral  in  Equation  (G-6)  is  trivial  and  is  given  by 


I.  - C2  [ (i*lAx  - iAx  ] + CD  [ ((i*l)Ax)2  - (iAx)2  ] 


♦ p-  [ ((i+1) Ax) 3 - (iAx) 3 ] 


CG- 19) 


The  complete  L2  norm  then  is 


l“-“NH0  ■ [ h * h * h l1* 


(G-20) 


L 


| 
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Appendix  H 

Computer-Generated  Plots  of  Results 

This  appendix  contains  the  graphical  results  of  the  project.  The 
results  are  handled  in  this  way  since  there  is  a large  number  of  plots 
which  need  to  be  carefully  organized  in  order  to  prevent  confusion. 

Each  three-letter  run  identifier,  for  example,  CER,  is  associated  with 
a data  set.  The  results  from  a given  set  of  data  may  be  displayed  in  a 
number  of  ways.  There  are  three  basic  formats  used.  The  first  shows 
the  discretization  error  ratio  as  a function  of  a . This  format 
emphasizes  the  effect  on  the  discretization  error  ratio  for  small  devi- 
ations in  a near  the  optimum  value.  It  also  shows  how  the  discretization 
error  ratio,  which  is  discussed  in  the  chapter  on  results,  varies  as 
a varies  over  its  range  of  stable  values.  The  results  from  three  time 
steps  are  given  in  each  case.  The  time  steps  shown  are  usually  the  last 
three  out  of  a total  of  eight  in  a given  calculation,  since  the  first  few 
time  steps  show  some  oscillatory  behavior,  but  then  smooth  out  in  later 
time  steps.  The  second  format  shows  the  discretization  error  ratio  as  a 
function  of  time.  Only  three  values  of  the  parameter  a were  followed 
in  time  using  this  format.  These  were  usually  the  values  for  the 
Crank-Nicolson,  optimum- imp licit,  and  pure-implicit  methods.  This 
second  format  shows  the  behavior  of  the  solution  for  the  predicted  values 
of  a , This  format  is  very  informative,  since  the  peak  in  a discreti- 
zation error  ratio  versus  a curve  may  appear  to  be  displaced  from  the 
optimum  value  of  a in  some  cases,  and  to  show  a high  convergence  rate. 
These  high  peaks  are  only  transient,  however,  and  are  superimposed  upon 


Li 
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the  true  peak,  whose  behavior  is  perhaps  better  indicated  by  the 
second  format.  The  third  format  shows  the  absolute  magnitude  of  the 
error  as  a function  of  ct  . The  discretization  error  ratio  is  a 


useful  tool  for  convergence  studies,  but  the  absolute  magnitude  of  the 
error  shows,  in  the  most  direct  way,  the  effect  of  the  parameter  a 
on  the  accuracy  of  the  solution. 

Each  of  these  three  formats  is  used  to  display  the  four  error 
norms  for  finite-elements  and  the  three  for  finite-differences.  Because 
of  the  large  amount  of  data  that  could  be  generated,  the  pointwise 
error  is  measured  only  at  the  point  x = 0.1  . The  generalized  mean 

is  the  sum  of  the  absolute  values  of  the  pointwise  errors  at  nine  (or  ten 
for  the  secondary  problem)  evenly  spaced  nodes.  This  shows  the  effect 
of  the  parameter  a on  the  pointwise  error  over  the  whole  interval.  The 
error  norm  termed  "maximum  error  any  node",  or  more  properly  called  a 
discrete  Tchebycheff  norm,  may  be  the  error  measure  of  greatest  interest 
to  the  engineer.  It  shows  the  effect  of  the  parameter  a on  the  maximum 
deviation  at  any  node  between  the  true  solution  and  the  numerical  solution 
The  norm  is  defined  only  for  the  finite-element  method  by 

llu-uNU0  3 t / (u-uV  dx  ]**  (H-l) 

o 

Since  the  numerical  solution  is  defined  for  all  points  on  the  interval  in 
the  finite-element  method,  not  just  at  the  nodes,  this  error  norm  shows 
how  the  accuracy  of  the  complete  numerical  solution  varies  with  the 
parameter  a . 
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In  the  analysis  of  the  primary  problem,  to  obtain  the  discretization 
error,  the  error  obtained  when  the  domain  was  divided  into  10  intervals 
was  compared  to  the  error  obtained  when  the  domain  was  further  subdivided 
into  20  intervals.  Then  the  error  for  a 20-interval  arrangement  was 
also  compared  to  the  error  for  a 40-interval  arrangement.  In  the  analysis 
of  the  secondary  problem,  only  the  10  to  20- interval  error  ratio  was 
computed.  The  more  complete  analysis  in  the  primary  problem  allows 
inspection  of  the  results  as  the  space  domain  is  further  subdivided.  The 
parameter  DX  in  the  legend  for  the  graphs  represents  Ax  . P is 
the  Fourier  modulus  A0/(Ax)2  and  a is  the  parameter  which  measures  the 
weight  placed  on  the  temperatures  in  the  new  time  step  in  the  numerical 
scheme,  or  in  other  words,  the  "degree  of  implicitness."  A value  of  0 
for  a would  correspond  to  Euler's  method,  or  for  finite-differences, 
an  explicit  scheme.  A value  of  .5  for  a gives  the  famous  Crank - 
Nicolson  scheme,  and  when  a is  1 , the  scheme  is  fully  implicit. 

In  the  graphs  for  the  finite-element  method,  results  are  shown 
only  for  values  of  a between  0.4  and  1.0  , since  for  the  values  of 
the  Fourier  modulus  used  in  the  investigation,  the  finite-element  method 
is  more  restrictive  with  respect  to  stability,  and  the  solutions  for  values 
of  a less  than  .5  are  often  subject  to  severe  oscillations  or  perhaps 
even  instability. 

Each  section  of  graphs  is  introduced  with  a descriptive  note  and  a 
key  to  aid  in  the  analysis.  In  each  section,  reference  will  be  made  to 
three  options  under  which  the  calculations  were  performed.  Option  one 
is  the  straightforward  approach  which  has  no  modifications.  In  option 
two,  the  exact  analytical  solution  has  been  substituted  at  the  first 
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tiae  step.  This  transforms  the  problem  into  a new  problem  without  the 
discontinuity  between  the  initial  and  boundary  conditions.  In  option 
three,  an  attempt  is  made  to  reduce  the  effect  of  the  discontinuity  in 
the  problem  by  subdividing  the  space-time  grid  for  the  first  time  step. 
The  most  extensive  analysis  is  made  with  option  one,  since  it  shows 
most  clearly  the  predicted  behavior  of  the  optimum  implicit  methods  for 
both  finite-difference  and  finite-elements. 
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Section  I 

The  Results  for  the  Primary  Problem  Using  Finite-Differences 

This  section  shows  the  graphical  results  for  the  solution  of 
the  primary  problem  by  finite-differences.  The  following  key  shows 
which  options  are  included  in  this  set  of  graphs. 


Table  H-I 

Key  to  the  Plots  in  Section  I 


Run  Identifier 

Fourier  Modulus  (p) 

Option  Number 

CDD 

0.5 

1 

1.0 

0 

1.0 

1 

CDG 

1.0 

2 

CDH 

2.0 

1 

OfT«t.  fTHISC  CfM.  I«l. 


FINITE 

DIFFERENCES 


<t.O  0.2 

0.4  0.6 

ALPHA 

--■&  - 

0.8 

1.0 

0X=  0.100 

TOi  0.050 

TIME 

0 

▲ 

0.030 

0.035 

P = 0.500 

LEGEND i 

+ 

0.040 

Fig.  H-l.  Discretisation  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  for  the  first  time  step. 
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DISC.  ERROR  RATIO 

cp.OO  2.70  5.40  8.10  10.79  13.49  16.19 

^ i ■ • ■ ■ a a 
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0X=  0.050 
TOt  0.025 
P = 0.500 


0.4  0.6 

ALPHA 

TIME 

LEGEND i 


0 0.030 
A 0.035 
+ 0.040 


Fig.  H-2.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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ERROR  RATIO 


DISC.  ERROR  RATIO 

cp.00  2.69  5.38  8.07  10.76  13.45  16.14 
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TIME 
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0 0.030 
A 0.035 
+ 0.040 


Fig.  H-4 . Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

TTie  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


142 


ERROR  RATIO 
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P = 0.500 

LEGEND i 
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scretization  Er 
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e step. 


*um*m 


Fig.  H-6.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 


solution  at  the  first  time  step. 
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CDD-OPT=l . PTH18E  ERR.  X=l. 


FINITE^-^ 

JW^ERENCES 


(P.OOO  0.010 

0.020 

TIME 

0.030 

0.040 

0X=  0.100 

RLPHfi 

q 0.3333 

TOt  0.050 

A 0.5000 

P = 0.500 

LEGEND t 

+ 1.0000 

Fig.  H-7.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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TWITE 
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1 

ft  000  0.010 

0.020 

0.030 
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OX-  0.050 

flLPHR 
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TOi  0.025 

A 0.5000 

P = 0.500 

LEGEND i 

+ 1.0000 

Fig.  H-8.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


ERROR  RATIO 


Fig,  H-9.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 


solution  at  the  first  time  step. 
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ALPHA 
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A 0.5000 

P = 

0.500 

LEGEND: 

+ l . 0000 

Fig.  H-10.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-ll . Discretization  Error  R 
The  exact  solution  has 
solution  at  the  first 


C0D-0PT=1.  OEM  HERN  ERROR 
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TOi  0.025 
P = 0.500 
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A 0.5000 
+ 1.0000 


0.040 


Fig.  H-12.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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COD-OPT=l.  PTW13E  ERH.  X=l. 
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Fig.  H-13.  The  Absolute  Magnitude  of  the  Discretization  Error  Versus 
Alpha  for  Problem  One.  The  exact  solution  has  been 
substituted  for  the  numerical  solution  at  the  first 
time  step. 
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ABSOLUTE  ERROR 


Fig.  H-  15.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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DISC.  ERROR  RATIO 

cP.OO  17.54  35.09  52.63  70.18  87.72  105.  2 

w a ■ i i i ■ I 
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0 0.060 
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Fig.  H-22.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


H-23.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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OISC.  ERROR  RATIO 

eP.OO  2.41  4.83  7.24  9.86  12.07  14.49 


Fig.  H-24 . Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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OISC.  ERROR  RATIO 

eP.OQ  2.69  5.37  8.06  10.75  13.44  16.12 

a i a a i i I 
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Fig.  H-25.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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F ig.  H-27.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


DISC.  ERROR  RATIO 

FD.OQ  2.42  4.84  7.27  8.69  12.11  14.53 


Fig.  H-28.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


DISC.  ERROR  RATIO 

53.00  2.70  5.39  8.09  10.78  13.48 


I C0f-0PT=1.  PTHISE  EHW>  X»t. 


TOt  0.025  A 0.5000 

P * 1.000  LEGEND*  + UQ000 


Fig.  H-29.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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DISC.  ERROR  RATIO 

50.00  2.41  4.83  7.24  9.66  12.07  14.49 
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Fig.  H-30.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-31.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-32.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


i» ion  | rror  Ratio  Versus  Time  for  Problem  One. 
v -.4.  t ». Motion  ha«  been  substituted  for  the  numerical 
m i. >11  •*■  at  the  first  time  step. 
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Fig.  H-34.  The  Absolute  Magnitude  of  the  IHscretization  lirror 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 


m 


Fig.  H-35 . The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-36.  The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-37.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-39.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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DISC.  ERROR  RATIO 

cP.OO  0.91  1.82  2.73  3.64  4.55  5.46 

w i i i i > i i 


cdo-opt=2.  nnx  error  any  node 


Fig.  H-40.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-41.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One 
The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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DISC.  ERROR  RATIO 
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Fig.  H-42.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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DISC.  ERROR  RATIO 
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Fig.  H-43.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-44.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-45.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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ERROR  RATIO 
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Fig.  H-46.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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AD-A056  508  AIR  FORCE  INST  OF  TECH  WR I GHT -PATTERSON  AFB  OHIO  SCH— ETC  F/S  20/13 
AN  INVESTIGATION  OF  THE  NUMERICAL  METHODS  OF  FINITE  DIFFERENCES — ETC(U) 
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DISC.  ERROR  RATIO 
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Fig.  H-47.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


185 


i COW-OPT»l.  OEM  r.gqH  EWROK 


ALPHA 

TIME  O 

A 0.140 
LEOENO.  + l60 

or  Ratio  Versus  Alpha  for  Problem  One. 
i has  been  substituted  for  the  numerical 
rst  time  step. 


186 


DISC.  ERROR  RATIO 

5D.00  3.96  6.72  10.09  13.45  16.81  20.17 


COH-Of T=1 . fTMiaC  X=. I 


FINITE 

DIFFERENCES 


0.040 


0X=  0.100 
TO*  0.050 
P = 2.000 


0.060  0.120 

TIME 

ALPHA  © 0.4583 
A 0.5000 
LEGENO:  + i.QOOO 


0. 160 


Fig.  H-49.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-50.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-Sl.  Discretization  F.rror  Ra 
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solution  at  the  first  t 
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53.00  2.70  5.41  8.11  10.82  13.52  16.29 


COH-OPT*!.  HA*  ERROR  ANY  NOOE 


"Finite 

DIFFERENCES 


000 


0.040 


ox-  0.050 
TOt  0.025 
P = 2.000 


0.080 

TIME 


ALPHA 


LEGEND i 


0.120 

q 0.4583 
A 0.5000 
+ 1.0000 


0.160 


Fig.  H-52.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-54.  Discretization  F.rror  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-55.  The  Absolute  Magnitude  of  the  Discretisation  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
< first  time  step. 
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Fig.  H-5b.  The  Absolute  Magnitude  of  the  Discretisation  F.rror  Versus 
Alpha  for  Problem  One.  The  exact  solution  has  been 
substituted  for  the  numerical  solution  at  the  first  time 
step. 


Fig.  H-57 . The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Section  II 

The  Results  for  the  Primary  Problem  Using  Finite-Elements 

This  section  shows  the  graphical  results  for  the  solution  of 
the  primary  problem  by  finite-elements.  The  following  key  shows 
which  options  are  included  in  this  set  of  graphs. 


Table  H-II 

Key  to  the  Plots  in  Section  II 


Run  Identifier 

Fourier  Modulus  (p) 

Option  Number 

CER 

0.5 

1 

CES 

1.0 

0 

CET 

1.0 

1 

CEU 

1.0 

2 

CEV 

1.0 

3* 

CEY 

2.0 

1 

* This  non  was  made  using  the  modification  used  in  Figs.  13  and  14 
in  Chapter  III.  This  modification  was  an  attempt  to  better 
approximate  the  initial  conditions. 
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Fig.  H-58.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-59.  Discretization  F.rror  Ratio  Versus  Alpha  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-60.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-61.  Discretization  Error  Ration  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-62.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-63.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-66.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-67 . Discretization  Error  Ratio  Versus  Time  for  Pi'oblem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-68.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 


solution  at  the  first  time  step. 
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Fig.  H-69.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 


solution  at  the  first  time  step. 
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Fig.  H-70.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-71.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-73.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-74.  The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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H-75 . The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Fig.  H-76.  The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Fig.  H-77. 


The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-79.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 
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H-8S.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 
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Fig.  H-86.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-88.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-89.  Discretization  lirror  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Pig.  H-90.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-91.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The. exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-92.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-93.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 


The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-94 . Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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5®.  00  2.70  5.39  8.09  10.78  13.48  16.17 


CET-OPT*l.  PTWUE  EM.  X-.1 


riNEAR  FINITE 
ELEMENTS 


0.020 


0X=  C.050 
TOi  C.025 
f * 1.000 


0.040 

TIME 


RLPHR 


LEGEND i 


0.060 

0 0.5000 
A 0.5833 
+ 1.0000 


0.080 


Fig.  H-95 . Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-96.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-97.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-98.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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DISC.  ERROR  RATIO 
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Fig.  H-102.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-103.  The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 


Fig.  H-104.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-105.  The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Fig.  H-106.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 

time  step  to  obtain  a more  accurate  solution.  ^ 
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Fig.  H-107.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first, 
time  step  to  obtain  a more  accurate  solution. 


Fig.  H-108 . Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-109.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

TTie  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-110.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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DISC.  ERROR  RATIO 

cp.OO  0.86  1.72  2.58  3.44  4.31  5.17 


Fig.  H-lll.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-112.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
tiae  to  obtain  a more  accurate  solution. 
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DISC.  ERROR  RATIO 
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Fig.  H-113.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-114.  Discretization  F.rror  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  have  been  modified  by  the  technique 


used  to  generate  Figs.  15  and  16. 
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DISC.  ERROR  RATIO 

cP.OO  0.78  1.57  2.35  3.14  3.92  4.70 


cp.00  0.80 


Fig.  H-116.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  have  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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OISC.  ERROR  RATIO 

cp.OO  0.77  1.55  2.32  3.10  3.87  4.65 


Fig.  H -117.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  ahve  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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DISC.  ERROR  RATIO 

JJ.00  0.79  1.58  2.38  3.17  3.96  4.75 


CEV-0PT=3.  OEN  MEAN  ERROR 


Fig.  H-118.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  have  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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DISC.  ERROR  RATIO 

cp.00  0.77  1.54  2.32  3.09  3.86  4.63 


Fig.  H-119.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  have  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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Fig.  H-120.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  conditions  have  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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Fig.  H-121.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  initial  condition  has  been  modified  by  the  technique 
used  to  generate  Figs.  15  and  16. 
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ERROR  RATIO 


Fig.  H-122.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-123.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-124.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-125.  Discretization  Error  I 
The  exact  solution  ha: 
solution  at  the  first 
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Fig.  H-126.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-127.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Ftg  H-IJ*.  Discretization  F.rror  Ratio  Versus  Alpha  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


OISC.  ERROR  RATIO 

cp.00  0.87  1.74  2.60  3.47  4.34  S.2L 


Fig.  H-130.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-131.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-132.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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DISC.  ERROR  RATIO 

50.00  2.70  5.41  8.11  10.82  13.52  16.23 


Fig.  H-133.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-134.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-136.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-137.  Discretization  Error  Ratio  Versus  Time  for  Problem  One. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig-  H-138.  The  Absolute  Magnitude  of  the  Discretization  Error 


Versus  Alpha  for  Problem  One.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Fig.  H-140.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Troblem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 


Fig.  H-141.  Hie  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  One.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Section  III 

The  Results  for  the  Secondary  Problem  Using  Finite-Differences 

This  section  shows  the  graphical  results  for  the  solution  of 
the  secondary  problem  by  the  method  of  finite-differences.  The 
following  key  shows  which  options  are  included  in  this  set  of  graphs. 


Table  H-III 

Key  to  the  Plots  in  Section  III 


DISC.  ERROR  RATIO 

cP.OO  4.66  9.33  13.99  18.65  23.31  27.98 
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Fig.  H-145.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-146.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-147.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-148.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-149.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-150.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-151. 


The  Absolute  Magnitude  of  the  Discretization  Error 
Versus  Alpha  for  Problem  Two.  The  exact  solution 
has  been  substituted  for  the  numerical  solution  at 
the  first  time  step. 
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Fig.  H-152.  The  Absolute  Magnitude  of  the  Discretisation  Error  Versus 
Alpha  for  Problem  T\vo.  The  exact  solution  has  been 
substituted  for  the  numerical  solution  at  the  first  time 
step. 
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Fig.  H-153.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  Two.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 


Fig.  H-154.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 


The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 

293 


OISC.  ERROR  RATIO 

J.  00  1.50  2.99  4.49  5.99  7.49  8.98 


COO- Of let.  MX  ERROR  ANT  NODE 


Fig.  H-155.  Discretization  Error  Ratio  Versus  Alpha  f 
The  space-time  grill  has  been  subdivided  f 
time  step  to  obtain  a more  accurate  solut 
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Fig.  H -156.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 


295 


Section  IV 

The  Results  for  the  Secondary  Problem  Using  Finite-Elements 

This  section  shows  the  graphical  results  for  the  solution  of  the 
secondary  problem  by  the  method  of  finite-elements.  The  following  key 
shows  which  options  are  included  in  this  set  of  graphs. 


Table  H-IV 

Key  to  the  Plots  in  Section  IV 


Run  Identifier 

Fourier  Modulus  (p) 

Option  Number 

CES 

1.0 
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1.0 
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CEU 

1.0 

2 
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cea-opr=o.  pTHiae  err.  x=. i 


Fig.  H-  157.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 
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Fig.  H>160.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  'IVo. 


Fif.  H-161.  Discretisation  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-162.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-163.  Discretiiation  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-164.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 


The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-165.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 


The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-166.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H-167.  Discretization  Krror  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 


DISC.  ERROR  RATIO 

50.00  0.87  1.74  2.60  3.47  4.34  5.21 


Fig.  H-168.  Discretization  Error  Ratio  Versus  Time  for  Problem  Two. 

The  exact  solution  has  been  substituted  for  the  numerical 
solution  at  the  first  time  step. 
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Fig.  H -169.  The  Absolute  Magnitude  of  the  Discretization  Error  Versus 


Alpha  for  Problem  TVo.  The  exact  solution  has  been 
substituted  for  the  numerical  solution  at  the  first 
time  step. 
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Fig.  H-170.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  Two.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 
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Fig.  H-171.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  Two.  TTie  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 


first  time  step. 


ROR 


Fig.  H-172.  The  Absolute  Magnitude  of  the  Discretization  Error 

Versus  Alpha  for  Problem  Two.  The  exact  solution  has 
been  substituted  for  the  numerical  solution  at  the 
first  time  step. 

312 


ERROR  RATIO 


CEU-OM«t.  RTMISC  ERR.  X«.t 


LINEAR  FINITE 
ELEMENTS 


**.4 

0.5  0.6 

0.7  0. 

8 

0.9 

1.0 

ALPHA 

OX-  0.100 

TINE 

0 

0.060 

TOt  0.050 

▲ 

0.070 

P s 1.000 

LEGEND i 

+ 

0.080 

Fig.  H-173.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-174.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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Fig.  H-176.  Discretization  Error  Ratio  Versus  Alpha  for  Problem  Two. 

The  space-time  grid  has  been  subdivided  for  the  first 
time  step  to  obtain  a more  accurate  solution. 
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The  transient  heat  conduction  equation,  with  Dirichlel  and  Neumann  boundary 
conditions,  is  solved  by  the  methods  of  finite-differences  and  finite-el ere 
and  the  numerical  solutions  are  investigated  with  respect,  to  accuracy  and 
stability.  A general  six  point  finite-difference  expression  is  used'  for 
which  there  exists  a high  order  a c curs. to  modification.  The  finite- element 
method  used  is  based  on  a stationary  variational  principle,  for.  -nl  method 
for  treating  accuracy  and  convergence  problem-  which  result  from  a discof  - 
5 tv  in  the  initial  condition  are  j nv-'s  t i •a"  • 1.  The  -.’rank-"  1 eel  cot.  ■ *1 
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is  n special  case  of  both  the  f inite-differor.ce  and  finite-element  " • the  . 

: • n fferenc®  version  of  the  Crank-1  ioolaon  b hod 
rare  accurate  thin  the  finite-element  version,  especially  when  a discon4  inui t; 

exists  1 x ti  ' Ltial  - q 1 : ti  on  an  I b bou >onditiona.  Dm  ' jh 

order  accurate  sc!  -<ros  for  both  finite— difference;!  arid  f i ni t e-elements  are 
g|  vn  * c i ■ ■ tuival  a I for  t 1 ,o  case  of  i ■ c 1 ‘ ■ . Some  oi  the  resu  s 

suggest  the  possibility  of  finding  a finite-element,  scheme  which  is  highly 
accurate  ■ r.»-ar  s nai'*e  sense  over  ! c at  ire  a Lution  doma in* 
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